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1  Proposed  Research  and  Technical  Objectives 
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instantaneous  aerodynamic  forces  and  therefore  flight  behavior.  The  wing  shape  and  how  it 
moves  through  the  flow  is  the  primary  mechanism  for  momentum  transfer  to  the  fluid  in  the 
flapping  flights.  During  this  project  period,  high-fidelity  CFD  simulation  and  adjoint 
analysis/optimization  tool  has  been  used  to  study  the  aerodynamic  role  of  wing  morphology 
models  built  from  the  previous  year  and  their  interaction  in  a  vast  parametric  space  including 
aspect  ratio,  direction  of  flexibility  (chord-wise,  span-wise,  etc.).  Aerodynamic  performance  has 
been  evaluated  by  looking  at  lift  production,  power  consumption,  and  efficiency.  Novel  tools  for 
studying  wing  morphing  during  complicated  flapping  flights  have  been  developed  to  understand 
the  underlying  physics  of  flexible  wings  in  flying  insects  and  birds  towards  the  bio-inspired  wing 
designs  with  superior  aerodynamic  performance.  In  this  funded  work,  we  built  up  a  team  with 
experts  from  two  institutes  and  hope  to  address  the  following  basic  questions  about  the 
aerodynamics  of  flapping  flight: 

1.  What  are  the  detailed  features  of  morphing  flapping  wings  and  how  to  build  up  high-fidelity 
kinematic  models? 

2.  What  is  the  proper  way(s)  to  characterize  the  wing  morphing  and  describe  the  key  motion 
components? 

3.  What  is  the  vortex  dynamics  and  associated  aerodynamics  in  morphing  wing  flapping  flight? 

4.  How  to  design  morphing  wing  kinematics  to  achieve  the  optimal  performance  of  the  wing 
beyond  biology? 

5.  What  is  optimal  solution  for  2D  and  3D  rigid/flexible  flapping  wings,  and  what  we  can  learn 
from  these  mathematically  optimal  solutions? 


2  Physics-based  Morphology  Analysis 


In  this  part,  four  sections  regarding  the  physics-based  morphology  analysis  on  the  proposed 
topic  are  presented.  We  briefly  summarize  these  works  in  the  following. 

(1)  Joint-based  morphing  surface  reconstruction 
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>  A  high-accurate  kinematics  measurement  tool  is  developed  to  quantify  the  propulsor 
flexion  and  body  kinematics  of  animals  in  free  flight. 

(2)  High-fidelity  Numerical  Simulations  on  the  Aerodynamics  of  Morphing  wings 

>  An  immersed  boundary  method  for  deformable  attaching  bodies  (IBM-DAB)  is 
developed  to  handle  direct  numerical  simulations  in  some  extreme  situations  that  are 
commonly  exist  in  nature,  including  solid  bodies  with  sharp  edges  and  those  with 
deformable,  attaching  membrane  bodies. 

>  The  3D  wake  structures  and  aerodynamic  performance  of  a  freely  maneuvering 
hummingbird  is  studied  in  detail.  Our  simulation  results  show  asymmetric  wake  structure 
between  inner  and  outer  wings  of  the  hummingbird.  A  unique  duel-ring  vortex  structure, 
which  is  the  source  of  the  wake  asymmetry,  is  found  in  the  wake  of  one  of  the  two  wings 
of  the  hummingbird.  The  duel-ring  vortex  structure  corresponds  to  larger  wing  twisting 
and  lower  drag  production,  which  creates  unbalanced  aerodynamic  forces  to  help  with  the 
maneuver 

>  Free  forward  flight  of  cicadas  has  been  investigated  through  high  speed  photogrammetry, 
three-dimensional  surface  reconstruction  and  computational  fluid  dynamics  simulations. 
We  report  two  new  vortices  generated  by  the  cicada  body.  One  is  the  thorax  generated 
vortex,  which  helps  the  downwash  flow,  indicating  a  new  phenomenon  of  lift 
enhancement.  Another  is  the  cicada  posterior  body  vortex,  which  entangles  with  the 
vortex  ring  composed  of  wing  tip,  trailing  edge,  and  wing  root  vortices. 

>  A  new  lift-enhancement  mechanism  due  to  the  wing-body  interaction  has  been  found  in 
cicada  forward  flight.  Our  results  have  shown  that  due  to  wing-body  interaction,  the 
wing-body  model  of  the  cicada  flight  had  a  18.7%  increase  in  total  lift  production 
compared  with  the  sum  of  the  lift  generated  in  body-only  and  wings-only  models,  and 
about  65%  of  this  enhancement  was  attributed  to  the  body.  This  resulted  from  a  dramatic 
improvement  of  body  lift  production  from  2%  to  11.6%  of  the  total  lift  produced  by  the 
wing-body  system.  Further  analysis  of  the  associated  near-field  and  far-field  vortex 
structures  has  shown  that  this  lift  enhancement  was  attributed  to  the  formation  of  two 
distinct  vortices  shed  from  the  thorax  and  the  posterior  of  the  insect,  respectively,  and 
their  interactions  with  the  flapping  wings. 

>  Aerodynamic  effects  of  the  wing  morphing  in  dragonfly  forward  flight  is  computationally 
investigated.  We  use  the  reconstructed  model  to  explore  the  effects  of  morphing  wings, 
first  by  removing  camber  while  keeping  the  same  time-varying  twist  distribution,  and 
second  by  removing  both  the  camber  and  the  spanwise  twist.  Our  simulation  results 
revealed  that  the  surface  deformation  can  improve  the  aerodynamic  functions  in  two 
ways:  1)  improving  the  power  economy  by  preventing  the  tip  vortex  bursting;  and  2) 
improving  the  leading-edge  vortex  attachment  by  suppressing  the  generation  of  the 
secondary  vortex.  As  a  result,  the  spanwise  twist  can  boost  the  aerodynamic  efficiency  up 
to  20%,  especially  during  the  wing  translational  phase. 

>  Aerodynamic  functions  and  associated  vortex  formation  of  the  two-pair  wings  of  a 
dragonfly  in  turning  maneuver  is  numerically  investigated.  The  wing  kinematics  analysis 
indicate  that  during  the  turn  there  are  large  asymmetries  between  the  wing  pitch  angle 
and  wing  stroke  angles  especially  for  forewings,  while  asymmetries  in  wing  deviation 
angle  are  generally  small.  The  asymmetrical  wing  kinematics  generate  unbalanced  forces 
in  both  vertical  and  horizontal  directions.  During  the  downstroke,  the  force  generated  by 
the  outer  wings  are  higher  than  inner  wings.  By  contrast,  the  inner  wings  generate  higher 
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forces  during  the  upstroke.  The  unbalanced  forces  lead  to  the  torques  along  its  body  axis 
in  three  directions.  Specifically,  both  forewings  and  hindwing  contribute  to  the  roll  torque 
generation,  whereas  forewings  are  dominant  in  generating  pitch  and  yaw  torques.  To 
examine  the  effect  of  forewing-hindwing  interaction  on  the  aerodynamic  performance, 
two  additional  cases,  one  with  forewings  only  and  the  other  with  hindwings  only,  are 
carried  out.  Our  results  indicate  that  the  interaction  can  lead  to  a  slight  forces  reduction, 
but  can  save  the  aerodynamic  power  up  to  11%. 

(3)  Low-dimensional  Modeling  (SVD)  and  Analysis  of  the  Morphing  wing  kinematics 

>  The  aerodynamic  roles  of  wing  morphing  of  a  hovering  dragonfly  is  studied  in  detail.  A 
spherical-coordinates-based  singular  value  decomposition  (SSVD)  method  is  developed 
and  applied  to  decompose  the  morphing  wing  kinematics  of  the  dragonfly  into  simple 
modes.  Results  have  shown  that  the  first  two  SSVD  modes  contain  93.1%  of  the  hovering 
wing  motion.  The  mode  1  (flapping  mode)  consists  of  a  simple  flapping  motion,  and  the 
mode  2  (morphing  mode)  contains  dynamic  wing  morphing  in  both  span-wise  and  chord- 
wise  directions.  By  evaluating  the  aerodynamic  role  of  the  SSVD  modes  using  a  high- 
fidelity  flow  simulation,  we  further  conclude  that  the  first  two  modes  can  recover  more 
than  96%  of  the  lift  production  and  91%  of  the  lift  economy  comparing  to  the  original 
flapping  wing  aerodynamics  whereas  the  mode  1  only  produces  5%  of  the  lift  and  4%  of 
the  lift  economy.  The  associated  flow  mechanisms  of  the  morphing  mode  are  found  to  be 
the  reduced  wing  tip  vortex  and  improved  attachment  of  leading-edge  vortex 

(4)  Optimization  of  low-dimensional  morphing  wing  models 

>  We  further  investigated  optimal  configurations  of  dominant  modes  on  aerodynamic 
performance  for  the  dragonfly  wing.  The  corresponding  optimized  low  dimensional  wing 
models,  which  can  beyond  biological  levels  of  aerodynamic  performance,  are  obtained. 
The  associated  flow  mechanisms  are  found  to  be  the  improved  LEV  strength  and  the 
reduced  TV  strength. 

2.1  Joint-based  3D  Surface  Reconstruction 

Recently,  advanced  photogrammetry  technology  has  been  used  to  study  flapping  wings  in 
nature.  Along  with  the  highly  accurate  surface  reconstruction  method  [1],  researchers  are  capable 
of  digitizing  detailed  propulsor  kinematics,  as  well  as  deformation  from  high-speed  images. 
However,  the  conventional  reconstruction  method  has  several  hard  constraints  regarding  some 
details  of  the  experiments,  such  as  the  camera  location/orientation,  the  lens  angle  of  view, 
marker  points  on  the  objects,  etc.  For  example,  it  is  usually  difficult  for  us  to  configure  three 
orthogonal  high-speed  cameras  to  conduct  experiments.  For  some  fast  and  agile  flyers,  such  as 
hummingbird,  it  is  not  possible  for  us  to  use  lenses  with  a  small  angle  of  view  due  to  their  large 
movements.  Perspective  error  becomes  a  big  issue  in  those  cases.  Moreover,  placing  marker 
points  on  the  wings  could  be  hard  due  to  the  wings’  size  or  properties.  A  more  robust 
reconstruction  method  is  needed  to  deal  with  diverse  flapping  wings  in  nature. 

A  joint-based  surface  reconstruction  method  is  developed  in  Autodesk  Maya  to  solve  the 
above  problems  in  this  work.  Here  we  use  a  case  of  a  maneuvering  hummingbird  to  demonstrate 
the  new  method. 

1 )  Static  Model  and  Virtual  Joints 
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In  order  to  perform  the  reconstruction,  a  realistic  static  model  of  the  hummingbird  with  some 
carefully  configured  motion  controllers,  which  are  so  called  “virtual  joints”,  is  needed  at  first. 


Figure  1.  (a)  Static  model  of  hummingbird  with  virtual  joints;  (b)  virtual  joints  configuration 
based  on  real  hummingbird  skeletal  anatomy. 


Figure  1(a)  shows  the  static  model  of  the  hummingbird  with  virtual  joints,  which  is  manually 
created  with  Autodesk  Maya.  The  dimensions  of  the  model  strictly  match  with  the  real 
hummingbird  we  are  going  to  study.  The  virtual  joints  have  six  degrees  of  freedoms  (three 
translations  and  three  rotations)  to  control  the  motion  of  adjacent  vertices  of  the  model.  Figure 
1(b)  shows  the  virtual  joints  configuration  at  the  hummingbird  wing,  which  presents  large 
deformation  during  the  flapping  motion.  This  kind  of  configuration  can  model  both  the  span  and 
chordwise  morphing  of  the  wing.  As  shown  in  the  figure,  the  virtual  joints  are  built  based  on  a 
real  hummingbird’s  skeletal  anatomy.  However,  some  parts  of  the  wing,  which  presents  large 
motion  in  flight,  do  not  really  have  bone  structures,  such  as  the  wing  tip  region.  We  also  add 
virtual  joints  to  control  the  motion  for  those  parts. 

2}  Wing  and  Body.  Surface  Reconstruction 


The  hummingbird  maneuvering  flight  is  then  reconstructed  using  the  joint-based  hierarchical 
subdivision  surface  method.  The  hummingbird  wings  are  marked  with  white  marker  points  to 
facility  the  three-dimensional  surface  reconstruction.  After  the  videotaping  was  done,  the  pose  of 
the  model  is  adjusted  to  match  one  frame  of  the  three  directions  of  high-speed  videos  by 
controlling  the  virtual  joints  in  six  degrees  of  freedom,  including  three  rotations  and  three 
translations.  Marker  points  on  the  wings  served  as  references  to  further  tune  the  location  of 
vertices  on  the  wing  models.  The  wing  models  are  generated  with  Catmull-Clark  subdivision 
surfaces  [2],  which  is  a  specific  cubic  spline  surface  representation  that  can  generate  smooth 
surfaces  from  meshes  of  arbitrary  topology  [3].  Similar  procedure  is  applied  to  the  hummingbird 
model  frame  by  frame.  As  an  example,  Figure  l(a,  b,  and  c)  show  the  poses  of  the  corresponding 
reconstructed  hummingbird  model. 
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(a)  (b)  (c) 


Figure  2.  Reconstructed  model  of  hummingbird  maneuvering  flight  at  (a)  t=0  ms,  (b)  t=47  ms, 
and  (c)  t=77  ms. 

2.2  High-fidelity  Numerical  Simulations  on  the  Aerodynamics  of  Morphing  Wings 

1 )  Immersed  Boundary  Method  for  Deformable  Attaching  Bodies 


The  governing  equations  considered  are  the  3D  unsteady  Navier-Stokes  equations  for  a 
viscous  incompressible  flow  with  constant  properties  given  by: 


dxi 

&U,  [  d 
dt  dxj 


1  dui  +  8  ^  du{  ' 

p  dt  dxj  ^  dxj 


(1) 


where  i,  j  - 1,2,3,  ut  are  the  velocity  components,  p  is  the  pressure,  and  p  and  v  are  the  fluid 
and  kinematic  viscosity. 

The  Navier-Stokes  equations(l)  are  discretized  using  a  cell-centered,  collocated  (non- 
staggered)  arrangement  of  the  primitive  variables  ui  and  p  .  In  addition  to  the  cell-center 


velocities  m,  ,  the  face-center  velocities  Ui ,  are  computed  (see  Figure  3).  The  equations  are 


integrated  in  time  using  the  fractional  step  method,  which  consists  of  three  sub-steps.  In  the  first 
sub-step  of  this  method,  a  modified  momentum  equation  is  solved  and  an  intermediate  velocity 
u  is  obtained.  A  second-order,  Adams-Bashforth  scheme  is  employed  for  the  convective  terms 
while  the  diffusion  terms  are  discretized  using  an  implicit  Crank-Nicolson  scheme,  which 
eliminates  the  viscous  stability  constraint.  In  this  sub-step,  the  following  modified  momentum 
equation  is  solved  at  the  cell-nodes: 
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Figure  3.  Schematic  describing  the  naming  convention  and  location  of  velocity  components 
employed  in  the  spatial  discretization  of  the  governing  equations. 


The  second  sub-step  requires  the  solution  of  the  pressure  correction  equation: 


n+ 1 
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At  p  Sxj 

which  is  solved  with  the  constraint  that  the  final  velocity  w"+1  be  divergence-free.  This  gives  the 
following  Poisson  equation  for  the  pressure  correction: 
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and  a  Neumann  boundary  condition  imposed  on  the  pressure  correction  at  all  boundaries. 

This  Poisson  equation  is  solved  with  a  highly  efficient  geometric  multigrid  method  [5], 
which  employs  a  modified  strongly  implicit  procedure  (MSIP)  [6]  smoother.  The  ability  to 
employ  such  methods  is  another  key  advantage  of  the  current  Cartesian  grid  approach  over  body- 
conformal  unstructured  grid  approaches.  Geometrical  multigrid  methods  are  relatively  simple  to 
implement  and  have  very  limited  memory  overhead.  Furthermore,  when  coupled  with  powerful 
smoothers  like  line-Gauss-Siedel  or  MSIP,  they  can  lead  to  a  numerical  solution  of  the  pressure 
Poisson  equation  which  scales  almost  linearly  with  the  number  of  grid  points.  In  contrast,  for 
unstructured  body-conformal  methods,  one  has  to  either  resort  to  algebraic  multi  grid  methods  [7, 
8]  or  other  more  complex  methods  such  as  agglomeration  multigrid  [9].  Another  choice  for 
solving  the  pressure  Poisson  equation  would  be  Krylov  subspace  based  methods  (such  as 
conjugate  gradient  or  GMRES)  but  these  require  effective  preconditioners  to  provide  good 
performance.  Our  past  experience  with  both  stationary  and  non-stationary  iterative  methods  [10- 
12]  indicates  that  geometric  multigrid  methods  are  very  well  suited  for  sharp  interface  immersed 
boundary  methods,  and  we  have  therefore  used  this  method  in  the  current  solver. 

Once  the  pressure  correction  is  obtained,  the  pressure  and  velocity  are  updated  as: 
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The  motions  of  natural  animals  are  reconstructed  using  the  method  introduced  in  Section  2.1.  As 
shown  in  Figure  4(a),  the  surface  of  a  reconstructed  model  is  represented  as  unstructured  mesh 
with  triangular  elements.  The  trunk  body  of  the  hummingbird  is  modeled  as  “solid  body”,  and 
the  wings  are  modeled  as  “membrane  body”,  which  presents  zero  thickness  in  the  simulation. 
The  two  types  of  bodies  firmly  attached  with  each  other,  and  the  boundary  conditions  near  the 
body  intersections  are  difficult  to  deal  with.  Here,  a  new  approach  is  introduced  to  overcome  the 
computational  difficulties. 
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Figure  4.  (a)  Example  of  a  hummingbird  model  with  unstructured  surface  mesh  of  triangular 
elements;  (b)  2D  schematic  describing  the  methodology  which  handles  deformable  attaching 
bodies. 


Figure  4(b)  shows  a  2D  schematic  plot  of  the  approach.  The  dashed  lines  are  immersed 
boundaries  and  the  solid  lines  are  the  corresponding  stair  boundaries.  For  the  membrane  body, 
the  stair  boundaries  are  chosen  as  the  closest  fluid  cell  faces  to  the  elements  on  the  immersed 
boundary.  As  shown  in  the  figure,  several  conflicting  faces  can  be  identified.  The  conflicting 
faces  are  stair  boundaries  which  directly  connect  to  the  intersection  of  the  solid  and  the 
membrane  stair  boundaries.  The  boundary  conditions  on  the  conflicting  faces  need  special 
treatments.  The  velocity  of  one  conflicting  face  may  be  calculated  according  to  a  reference 
velocity  of  the  solid  immersed  boundary  at  one  time  step.  For  the  next  time  step,  the  reference 
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velocity  may  be  altered  to  the  velocity  of  the  membrane  immersed  boundary  due  the  motion  of 
the  moving  boundaries.  This  will  cause  a  discontinuity  of  the  velocity  boundary  condition  on  the 
conflicting  face  for  those  two  time  steps,  which  will  cause  numerical  instability  and  eventually 
lead  to  the  divergence  of  the  flow  solver. 

In  order  to  solve  the  above  problem,  the  velocities  of  the  conflicting  faces  are  corrected 
accordingly.  Here  we  describe  the  method  using  the  upper  conflicting  face  (Figure  4b)  as  an 
example.  From  the  center  of  the  conflicting  face,  we  first  drop  two  lines  perpendicular  to  the 
solid  and  membrane  immersed  boundaries.  There  we  get  two  boundary  intercept,  BIs  and  BIm . 


Another  two  fluid  cell  faces  which  are  close  to  the  conflicting  face  can  be  found,  and  the 
corresponding  face  centers  are  Fj  and  F2 .  An  interpolation  stencil  (shaded  area)  is  formed 

according  to  the  four  points.  In  3D  cases,  the  stencil  consists  of  eight  points  (two  boundary 
intercept  and  six  fluid  cell  face  centers).  The  conflicting  face  velocity  can  then  be  calculated  as 
follows: 


(6) 


UCF  (xj ,  x2 ,  x3 )  =  Cjx jX2x3  +  C2x jX2  +  C3x2x3 

+C4XjX3  +  C5Xj  +  C6x2  +  C7x3  +  C8  +  0(A2  ^ 

The  eight  unknown  interpolation  coefficients  can  be  determined  in  terms  of  the  face 
velocities  of  the  eight  points  of  the  interpolation  stencil: 


{CVC2,...,C8}T  =[Vll  {Uv U2,...,U8}r  (7) 

where  Ui  are  the  face  velocities  at  the  interpolation  stencil  points  and  [V]  is  the  Vandermonde 
matrix.  This  process  continues  until  the  velocities  of  all  conflicting  faces  are  corrected. 


2 )  Asymmetric  Aerodynamic  Forces  and  Wake  Structures  of_a  Maneuvering  hummingbird 


Hummingbirds  perform  turning  maneuvers  as  often  as  they  hover  or  cruise,  especially  when 
they  need  to  forage  from  one  location  to  another.  However,  to  date,  turning  flight  has  received 
little  attention  and  there  is  no  detailed  forces  and  three-dimensional  flow  structure  data  to 
achieve  a  quantitative  analysis  of  hummingbird  in  a  turning  maneuver.  To  fill  this  gap,  a  high¬ 
speed  photogrammetry  system  and  three-dimensional  joint-based  surface  reconstruction 
technology  are  used  to  reveal  hummingbird  wing  kinematics  and  deformations  during  a  free 
maneuvering  flight.  The  aerodynamic  performance  is  then  studied  using  an  in-house  immersed 
boundary  method  (IBM)  based  computational  fluid  dynamics  (CFD)  solver.  To  the  best  of  our 
knowledge,  this  is  the  first  study  on  the  unsteady  aerodynamics  of  hummingbird  in  maneuvering 
flight. 
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Figure  5.  (a)  Schematic  plot  showing  hummingbird  body  yaw  motion;  (b)  time  course  of  body  of 
the  maneuvering  hummingbird.  A  top  view  of  the  hummingbird  at  the  top  of  the  figure  indicates 
the  yaw  throughout  the  maneuver.  Yaw  (red),  pitch  (black),  and  roll  (blue)  histories  are  shown 
first,  followed  by  path  position  histories  in  horizontal  (green)  and  vertical  (black)  direction. 

To  better  interpret  the  characteristics  of  body  and  wing  kinematics,  we  divided  the  whole 
maneuvering  process  into  three  phases  according  to  the  body  yaw  angle.  Figure  5  shows  a 
schematic  plot  of  hummingbird  body  yaw  motion  and  the  time  course  of  body  kinematics.  We 
can  see  from  Figure  5(b)  that  there  are  in  total  six  strokes  of  this  flight.  The  hummingbird  is  first 
at  an  “accelerating  phase”  (1st  stroke),  in  which  the  bird  initiates  the  yaw  turn.  The  body  yaw 
angle  start  to  increase  while  the  other  two  Euler  angles  stay  unchanged  within  this  phase.  After 
that,  the  hummingbird  enters  a  “turning  phase”  for  the  following  three  strokes  (2nd  to  4th  stroke). 
The  body  yaw  angle  keep  increasing  and  shows  an  oscillating  profile,  which  indicates  an  active 
control  of  the  turning  for  the  hummingbird.  The  body  roll  angle  shows  a  little  bit  decrease  and 
the  body  pitch  angle  stays  unchanged.  The  last  two  strokes  of  the  flight  (5  th  to  6th  stroke)  are 
called  the  “recovering  phase”,  in  which  the  body  yaw  angle  stays  at  a  high  value  and  shows  a 
small  increase,  while  the  other  two  Euler  angles  stay  unchanged.  The  hummingbird  is  recovering 
from  a  turning  status  to  a  hovering  status  within  this  phase.  The  path  position  of  this  flight 
(Figure  5b)  shows  a  little  bit  descending  motion  in  vertical  direction  for  the  first  two  strokes,  and 
the  motion  in  horizontal  direction  is  limited.  We  have  observed  similar  phases  in  other 
maneuvering  high-speed  videos  we  shoot.  The  similarity  suggests  that  the  body  motion  result 
from  similar  aerodynamic  or  dynamic  mechanisms. 

In  order  to  determine  how  hummingbird  change  wing  motion  to  perform  the  turning,  we  also 
study  the  wing  kinematics  of  this  flight.  Figure  6(a)  shows  the  three  Euler  angles  defining  the 
wing  position  in  the  wing-root  coordinate  system  (X'Y'Z'),  in  which  the  X'-axis  is  parallel  with 
the  body  longitudinal  direction,  the  F'-axis  is  along  the  lateral  direction  and  the  Z'-axis  complies 
with  the  right-hand  rule.  The  mean  stroke  plane  connected  the  wing  root  and  wingtips  at  the  start 
and  end  of  the  downstroke.  The  stroke  position  angle  cj)(t)  defines  the  angular  position  of  the 
wing  in  the  mean  stroke  plane,  with  0°  aligning  with  the  negative  direction  of  the  F'-axis.  The 
deviation  angle  0(t)  is  the  angle  between  the  base-to-wingtip  line  and  the  mean  stroke  plane.  The 
pitch  angle  a(t)  is  defined  as  the  angle  of  the  wing  chord  with  respect  to  the  tangent  of  the  wing 
trajectory. 
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Figure  6.  (a)  Schematic  plot  showing  hummingbird  wing  Euler  angles  definition;  Time  course  of 
wing  stroke  (b),  wing  deviation  (c),  and  wing  pitch  angles  (d)  during  the  turning  phase.  The  wing 
Euler  angle  histories  over  the  3rd  stroke  are  shown  in  solid  curves.  Angle  differences  for  other 
strokes  in  the  turning  phase  with  respect  to  the  3rd  stroke  are  shown  as  error  bars.  Red  and  blue 
correspond  to  the  inner  and  outer  wings,  respectively.  Shaded  areas  stand  the  downstrokes  and 
unshaded  areas  stand  the  upstrokes. 


The  time  course  of  wing  Euler  angles  during  the  turning  phase  are  shown  in  Figure  6(b,  c,  d) 
according  to  the  above  definitions.  We  can  see  from  the  figure  that  asymmetries  of  the  wing 
kinematics  between  inner  and  outer  wings  can  be  identified  in  the  turning  phase,  especially  for 
the  time  course  of  wing  deviation  angle  and  pitch  angle.  As  shown  in  Figure  6(b),  the  wing 
stroke  angle  history  shows  one  peak  for  a  flapping  cycle,  which  locates  at  the  end  of  the 
downstroke.  The  amplitude  of  the  wing  stroke  angle  is  around  120°  for  both  inner  and  outer 
wings.  However,  the  inner  wing  stroke  angle  is  a  little  bit  smaller  (7%  smaller)  during  the 
downstroke  in  turning  phase,  while  it  is  almost  the  same  comparing  to  that  during  the  upstroke. 
For  the  wing  deviation  angle  history,  the  plot  shows  two  valleys  for  a  flapping  cycle,  which  is 
located  near  the  mid  of  downstroke  and  the  mid  of  upstroke,  respectively.  Much  greater 
differences  can  be  observed  for  the  wing  deviation  angle  history.  The  inner  wing  deviation  angle 
history  shows  larger  amplitude  and  smaller  value  in  both  downstroke  and  upstroke.  This  is 
because  of  the  outer  wing  stroke  plane  of  the  hummingbird  is  tilted  up  in  the  turning  phase,  and 
also,  the  figure  eight  motion  of  the  outer  wing  is  smaller  comparing  to  that  of  the  inner  wing.  For 
the  wing  pitching  angle  histories,  all  plots  show  one  valley  and  one  peak  near  the  mid 
downstroke  and  mid  upstroke,  respectively.  More  importantly,  the  wing  pitching  angle  histories 
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show  great  asymmetry  between  the  inner  and  outer  wings  in  the  turning  phase.  The  valley  value 
of  the  inner  wing  pitching  angle  is  much  smaller  than  that  of  the  outer  wing,  while  the  peak  value 
is  much  greater.  It  leads  to  the  average  angles  of  attack  during  both  the  downstroke  and  the 
upstroke  being  much  smaller  for  the  inner  wing. 

CFD  simulation  is  conducted  using  the  kinematics  of  the  wings  and  the  body  from  the  three- 
dimensional  surface  reconstruction.  The  time  course  of  the  averaged  aerodynamic  forces 
produced  in  the  turning  phase  (average  over  the  2nd~4th  stroke)  are  shown  in  Figure  7(a,  b).  We 
can  see  from  Figure  7(a)  that  the  lift  coefficients  for  both  inner  wing  and  outer  wing  show  two 
peaks  in  one  complete  stroke.  The  peaks  locate  at  around  50%  of  each  half  stroke,  and  it  is  worth 
noting  that  the  average  values  of  the  lift  coefficient  in  downstroke  are  much  higher  than  that  in 
upstroke.  This  conforms  to  the  previous  conclusion  made  on  a  hovering  hummingbird  by  Song  et 
al.  [13].  However,  the  double  peaks  in  upstrokes,  which  are  also  reported  by  this  paper,  are  not 
observed  in  this  maneuvering  case.  A  possible  explanation  for  that  is,  due  to  the  body  motion  of 
the  maneuvering  hummingbird,  the  effects  of  wake  capture  during  the  stroke  reversal  is 
weakened. 


(a)  (b) 


Figure  7.  Time  course  of  lift  (a)  and  drag  (b)  coefficients  during  the  turning  phase  of 
hummingbird  pure  yaw  turn.  The  force  coefficient  histories  over  the  3rd  stroke  are  shown  in 
solid  curves.  Force  coefficient  differences  for  other  strokes  in  the  turning  phase  with  respect  to 
the  3rd  stroke  are  shown  as  error  bars.  Red  and  blue  correspond  to  the  inner  and  outer  wings, 
respectively.  Shaded  areas  stand  the  downstrokes  and  unshaded  areas  stand  the  upstrokes. 


Besides  the  lift  force,  horizontal  force  is  also  important  in  maneuvering  flight  since  it  can 
generate  torque  to  drive  the  turn,  especially  for  the  pure  yaw  turn  case.  Figure  7(b)  shows  the 
drag  coefficient  history  of  this  maneuvering  flight.  We  can  see  from  the  plot  that  the  horizontal 
force  generated  by  the  inner  wing  is  always  greater  than  that  generated  by  the  outer  wing.  Such 
force  asymmetry  in  horizontal  direction  can  accelerate  the  turn  in  downstroke  and  damps  the  turn 
in  upstrokes.  This  half-cycle-based  turning  strategy  is  more  flexible  comparing  to  some  other 
nature  flyers,  like  the  fruit  flies,  their  wings  always  generate  resultant  torque  towards  the  turning 
direction  when  they  perform  the  maneuvering  flight.  It  is  much  easier  for  the  flyer  to  stop  the 
turn  and  adopt  an  alternative  flight  motion  based  on  what  it  needs,  which  is  often  the  purpose  of 
pure  yaw  turn. 


17 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


The  vortex  formation  over  a  stroke  cycle  (the  3rd  stroke)  is  shown  in  Figure  8,  in  which  the 
vortex  structures  are  identified  by  the  Isosurface  of  Q-criterion  (Q=10).  The  vortex  structures  are 
colored  by  the  non-dimensional  pressure.  The  stroke  cycle  starts  from  the  downstroke.  Figure  8(a) 
shows  the  vortex  structures  at  t/T=0.24,  which  corresponds  to  the  peak  of  the  vertical  force 
during  downstroke  (Figure  8a).  The  most  significant  vortex  structures  at  this  moment  is  the 
formation  of  leading-edge-vortex  (LEV).  The  LEV  attaches  well  to  the  wing  surface.  More 
importantly,  the  LEV,  the  tip  vortex  (TV),  the  trailing-edge  vortex  (TEV),  and  the  root  vortex 
(RV)  are  connected  end  to  end,  forming  a  vortex  loop,  within  which  the  air  moves  downward.  A 
previously  shed  vortex  loop  (PVL)  can  also  be  observed.  It  was  generated  during  the  upstroke 
prior  to  this  moment,  and  the  PVL  is  connected  to  the  newly  formed  vortex  loop  by  the  TEV. 
Similar  vortex  structures  are  found  in  hovering  hummingbirds  in  previous  studies  [13-15]. 

As  time  advances  to  t/T=0.33,  as  shown  in  Figure  8(b),  which  corresponds  to  the  peak  of 
horizontal  force  in  downstroke  (Figure  8b),  the  wings  are  near  the  end  of  downstroke  and  rotate 
rapidly  along  their  own  axis.  The  outer  wing  LEV  is  divided  into  two  branches,  known  as  dual 
LEV  [3],  and  two  shed  vortices,  SLEV  and  STEV,  can  be  identified.  New  vortex  structures  can 
be  found  at  the  outer  side  of  the  hummingbird  for  this  time  instance.  The  LEV,  TV  and  STEV 
are  connected  with  each  other  to  form  a  vortex  ring  near  the  outer  wing  tip  region.  Also,  the  LEV, 
SLEV,  TEV  and  RV  are  connected  to  form  another  vortex  ring.  It  is  worth  noting  that  the 
directions  of  the  vortex  tube  STEV  and  SLEV  are  opposite.  Similarly,  the  PVL  can  be  observed 
and  it  is  connected  to  the  later  newly  formed  vortex  ring.  At  the  inner  side  of  the  hummingbird, 
the  LEV  develops  and  the  newly  formed  vortex  ring  grow  larger.  However,  the  major  vortex 
structures  stay  the  same  as  we  described  at  t/T=0.24. 

At  the  end  of  the  downstroke  (t/T=0.48),  as  shown  in  Figure  8(c),  the  dual-ring  vortex 
structure  at  the  outer  side  of  the  hummingbird  still  can  be  observed.  The  vortices  convect  further 
downstream  and  the  PVL  starts  to  dissipate.  The  dual-ring  vortex  structure  begins  to  merge  to  a 
single  vortex  ring. 

During  the  upstroke  (t/T=0.73),  as  shown  in  Figure  8(d),  which  corresponds  to  the  peak  of 
vertical  force  at  upstroke  (Figure  8a).  The  major  vortex  structures  are  similar  to  that  at  t/T=0.24, 
except  for  the  more  complex  and  stronger  vortices  due  to  larger  angle  of  attack  of  the  wings  in 
upstroke  (Figure  8c).  At  t/T=0.88,  which  corresponds  to  the  peak  of  horizontal  force  at  upstroke, 
the  dual-ring  vortex  structure  can  be  observed.  Differently,  it  is  the  inner  side  of  the 
hummingbird  that  presents  the  double  loop  vortex  structure.  At  the  end  of  upstroke,  as  shown  in 
Figure  8(f),  the  dual-ring  vortex  structure  still  exists  and  the  two  vortex  rings  start  to  merge  to  a 
single  vortex  ring. 

For  all  time  instances  discussed  above,  the  inner  side  and  outer  side  of  the  hummingbird 
show  significant  asymmetry  in  vortex  wake  structures.  The  unique  dual-ring  vortex  structure 
exists  in  the  outer  side  during  downstroke  and  exists  in  the  inner  side  during  upstroke.  Such  dual¬ 
ring  vortex  structure  is  responsible  for  the  wake  asymmetry.  The  kinematic  difference  of  the 
inner  and  outer  wings  may  result  in  the  unique  dual-ring  vortex  structure,  and  further  influence 
the  wake  structures  of  the  inner  and  outer  side  of  the  hummingbird. 
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Figure  8.  (a-f)  Time  course  of  vortex  development  of  the  3rd  stroke  of  the  hummingbird  pure 
yaw  turn,  visualized  by  the  Q-criterion.  The  vortex  structures  are  colored  by  non-dimensional 
pressure. 
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Concluding  Remarks: 

A  three-dimensional  direct  numerical  simulation  was  performed  for  a  maneuvering 
hummingbird  with  a  very  accurate  wing  and  body  motion  reconstructed  from  high-speed  images. 
The  simulation  captures  the  aerodynamic  force  and  power  characteristics  in  the  entire 
maneuvering  process  and  also  details  of  the  flow  field,  including  the  unsteady  wake  structures  in 
both  near  and  far  field.  Our  results  on  the  maneuvering  hummingbird  kinematics  suggest  that 
hummingbird  sustained  yaw  turns  by  tiling  up  the  outer  wing  stroke  plane  and  suppressing  the 
outer  wing  figure  eight  motion;  also  by  increase  the  wing  pitching  angle  in  downstroke  and 
decrease  the  wing  pitching  angle  in  upstroke  to  create  wing  angle  of  attack  asymmetry  between 
inner  and  outer  wings.  For  the  wake  structures,  strong  asymmetric  wake  topology  is  identified. 
More  importantly,  a  unique  dual-ring  vortex  structure,  which  is  the  source  of  the  wake 
asymmetry,  is  found  in  the  wake  of  one  of  the  two  wings  of  the  hummingbird.  The  dual-ring 
vortex  structure  corresponds  to  larger  wing  twisting  and  lower  drag  production,  which  creates 
unbalanced  aerodynamic  forces  to  help  with  the  maneuver. 

(3)  Computational  investigation  of  cicada  aerodynamics  in  forward  flight 

Many  aerodynamic  mechanisms  of  force  generation  by  flapping  wings  have  been  proposed 
based  upon  studies  on  rigid  mechanical  models,  including  wing-wake  interactions  and  rotational 
circulation [16],  delayed  stall  during  the  translation  portion  of  the  stroke  [17],  axial  flow 
stabilized  leading  edge  vortex  (LEV)  [18],  and  rotational  accelerations  [19].  Although  studies  on 
wing  models  have  substantially  advanced  the  understanding  of  insect’s  flight,  wing  model  has 
intrinsic  restrictions  such  as  wing  rigidity  and  simplified  kinematics.  In  this  work,  a  realistic 
wing-body  model  a  cicada  ( Tibicen  linnei )  is  employed  to  investigate  the  aerodynamic 
performance  and  the  wake  patterns  in  insect  forward  flight. 

The  cicada  flight  is  reconstructed  using  template-based  hierarchical  subdivision  surface 
method.  Marker  points  on  the  wings  are  digitized  in  each  recorded  image  at  each  time  step.  As 
an  example,  Figure  9a  shows  a  raw  picture  a  cicada  during  downstroke.  Figure  9b  presents  a 
comparison  between  a  real  cicada  and  the  reconstructed  result.  Figure  10  shows  the  three  Euler 
angle  quantifying  the  flapping  kinematics  of  the  forewing  and  hindwing. 


(b) 


(a) 


Leading  edee 


Figure  9.  Real  cicada  and  its  reconstruction,  (a)  Raw  picture;  (b)  Comparison  of  real  and 
reconstructed  cicada. 
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Figure  10.  Stoke  angle  (0°  lateral,  downstroke  positive),  deviation  angle  (upward  positive),  and 
pitch  angle  (rotation  around  wing  span,  smaller  than  90°  when  leading  edge  is  forward)  for  the 
forewing  (a)  and  the  hindwing  (b). 

After  the  model  is  built,  we  carry  out  the  high-fidelity  simulation  on  the  cicada  forward  flight. 
The  vortex  formation  over  a  stroke  cycle  is  shown  in  Figure  11,  in  which  the  vortex  structures 
are  identified  by  the  iso-surface  of  (9-criterion  ((9=10).  At  the  onset  of  the  downstroke  (Figure 
11a),  an  LEV  is  uniformly  distributed  along  each  leading  edge.  A  forewing  starting  vortex  (FSV) 
is  formed  along  the  extruded  forewing  trailing  edge  (partial  of  which  is  hooked  with  the  leading 
edge  of  the  hindwing).  The  LEV  connects  to  the  FSV  through  forewing  tip  vortex  (FTV).  On 
each  hindwing,  a  trailing  edge  stopping  vortex  (TSV)  is  created  and  is  about  to  be  shed  during 
the  pronation.  As  wings  further  stroke  down  (Figure  lib),  three-dimensional  LEV  develops  due 
to  the  onset  of  wingtip  vortices  [20],  A  hindwing  starting  vortex  (HSV)  forms  at  the  hindwing 
trailing  edge.  Surrounding  the  cicada  wing,  the  leading  edge  vortex,  forewing  tip  vortex,  and 
trailing  edge  vortices  (FSV  and  HSV)  all  together  form  a  horseshoe-like  vortex  structure.  At  the 
mid-downstroke  (Figure  lie),  the  LEV  is  fully  developed  and  a  cone  shape  FTV  is  formed  on 
the  dorsal  surface  of  the  wing.  The  starting  vortices  (FSV  and  HSV)  from  the  forewing  and  the 
hindwing  are  merged.  Near  the  root  of  the  hindwing,  a  wing  root  vortex  (WRV)  is  generated  on 
each  side  of  the  cicada.  At  the  end  of  the  downstroke  (Figure  1  Id),  the  vortex  rings  are  elongated. 
A  strong  tip  vortex  is  developed  from  the  hindwing  (HTV),  entangling  with  fore  wing  tip  vortices 
(FTV).  At  the  mid-upstroke  (Figure  lie),  the  vortices  emanated  from  the  posterior  cicada  body 
are  elongated  and  intertwined  with  the  wing  root  vortices.  A  stopping  vortex  (SV)  shed 
forwardly  during  supination  can  also  be  seen.  At  the  end  of  upstroke  (Figure  Ilf),  the  vortex  tube 
previously  formed  in  the  downstroke  from  the  posterior  cicada  body  has  been  detached,  with 
each  arm  hooked  to  the  vortex  ring  generated  in  the  downstroke.  Eventually,  the  vortex  ring  in 
each  side  of  the  body  will  merge  with  the  PBV,  and  potentially  form  one  big  vortex  ring  after 
several  stroke  cycles. 
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Figure  11.  Time  course  of  vortices  development,  visualized  by  Q-criterion.  The  left  and  middle 
columns  are  back  view  and  top  view  respectively,  and  the  vortex  structures  are  colored  by 
spanwise  vorticity.  The  right  column  shows  the  structures  in  projection  view,  colored  by 
streamwise  vorticity.  Color  bars  in  the  first  row  apply  to  figures  in  the  same  column. 


Figure  12a  shows  the  velocity  vector  and  two-dimensional  streamlines  in  a  transverse  plane 
cutting  through  the  wing  and  part  of  the  body.  The  LEV,  TXV  and  WRV  actually  form  a  jet 
through  which  strong  downwash  flow  is  induced  near  the  mid-wing  region.  Figure  12b  is  the 
velocity  vector  in  a  transverse  plane  half  the  body  length  behind  the  cicada,  at  which  the  trailing 
edge  stopping  vortices  (TS  V)  are  cut.  The  pair  of  TSV  approaches  to  the  center  plane  due  to  self- 
induction.  The  induced  flow  by  TSVs  can  potentially  enhance  the  formation  and  development  of 
the  wing  root  vortex,  and  therefore  indirectly  help  the  downwash.  The  wing  root  vortex  induces 
upwash  flow  proximal  to  the  hindwing  and  the  body,  but  the  affected  region  is  far  smaller  than 
the  estimated  upwash  region  generated  by  bumblebees  [21]. 


(a)  (b) 


Figure  12.  Transverse  plane  cut  at  mid-do wnstroke.  (a)  Cut  through  wing  and  body  (b)  Cut 
through  the  near  wake  (no  wings  or  body  being  cut).  Contour  of  Q-criterion,  velocity  vector,  and 
two-dimensional  streamline  seen  from  the  back  view. 


Concluding  Remarks: 

Free  forward  flight  of  cicadas  is  investigated  through  high  speed  photogrammetry,  three- 
dimensional  surface  reconstruction  and  computational  fluid  dynamics  simulations.  We  report 
two  new  vortices  generated  by  the  cicada  body.  One  is  the  thorax  generated  vortex,  which  helps 
the  downwash  flow,  indicating  a  new  phenomenon  of  lift  enhancement.  Another  is  the  cicada 
posterior  body  vortex,  which  entangles  with  the  vortex  ring  composed  of  wing  tip,  trailing  edge, 
and  wing  root  vortices.  Some  other  vortex  features  include:  independently  developed  left  and 
right  hand  side  leading  edge  vortex,  dual  core  leading  edge  vortex  structure  at  the  mid-wing 
region,  and  near  wake  two-vortex-ring  structure.  In  the  cicada  forward  flight,  approximately  79% 
of  the  total  lift  is  generated  during  the  downstroke.  Cicada  wings  experience  drag  in  the 
downstroke,  and  generate  thrust  during  the  upstroke.  Energetics  study  shows  that  the  cicada  in 
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free  forward  flight  consumes  much  more  power  in  the  downstroke  than  in  the  upstroke,  to 
provide  enough  lift  to  support  the  weight  and  to  overcome  drag  to  move  forward. 

(4)  New  lift  enhancement  mechanism  of  wing-body  interaction  in  cicada  forward  flight 

Following  the  above  study,  we  further  investigate  the  aerodynamic  role  of  those  body¬ 
generated  vortices  in  cicada  forward  flight.  A  new  lift  enhancement  mechanism,  which  is 
associated  with  wing-body  interaction,  is  found.  In  this  work,  in  order  to  examine  the  effects  of 
WBI,  two  simplified  models,  WN  (wings-only)  and  BD  (body-only),  are  created  based  on  the 
original  reconstructed  WB  (wing-body)  model.  In  the  WN  model,  the  same  wing  kinematics  are 
kept  as  the  WB  model.  For  the  BD  model,  only  the  body  with  the  same  inclination  angle  as  that 
of  the  WB  model  is  employed.  By  isolating  the  coupling  of  the  body  and  wing  models,  we  desire 
to  investigate  the  inherent  nature  of  body-involved  unsteady  force  generation  mechanism  by 
comparing  both  aerodynamic  performance  and  associated  wake  structures. 

The  comparisons  of  the  instantaneous  forces  on  the  cicada  body  and  the  right  wing  are  shown 
in  Figure  13(a)  and  (b),  respectively.  In  each  plot,  solid  lines  indicate  simulation  results  from  the 
WB  case  while  dashed  lines  represent  either  the  BD  or  WN  model.  The  overall  lift  force 
produced  by  the  wings  and  the  body  together  is  increased  by  about  18.7%  according  to  Table  1, 
which  suggests  a  significant  aerodynamic  benefit  generated  by  WBI.  This  is  mainly  because  of 
the  lift  increment  on  the  cicada  body,  which  contributes  about  65%  of  the  total  lift  enhancement. 
The  body  lift  accounts  for  about  1 1 .6%  of  the  total  lift  generation  in  the  WB  model,  whereas  the 
body  can  only  produce  less  than  2%  of  the  total  lift  when  separated  from  the  wings  and 
simulated  under  the  same  flow  conditions.  Therefore,  the  body  plays  a  more  important  role  in 
force  generation  when  it  interacts  with  the  flapping  wings. 


Table  1.  Cycle-averaged  lift  coefficient  and  its  enhancement  due  to  wing-body  interaction. 
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WB 
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Figure  13.  Time  traces  of  lift  and  thrust  coefficients  of  the  body  (a)  and  a  single  wing  (b)  during 
the  fourth  flapping  cycle  when  both  the  aerodynamic  forces  and  flow  reach  a  periodic  state.  Peak 
1  and  2  in  figure  (a)  occur  at  t/T=0.3  and  0.74,  respectively. 
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Pressure  distribution  on  the  body  and  the  wings  are  analysed.  The  pressure  coefficient  is 
defined  as  Cp  =  (p~ p,)/ 0.5 pU2 ,  where  pr  is  the  pressure  in  the  freestream.  The  above  force 

comparison  indicates  that  the  lift  on  the  body  increased  a  lot,  so  we  plot  the  surface  pressure  on 
the  body  in  Figure  14.  Figure  14a  is  the  pressure  on  the  body  in  the  WB  model  at  t/T=0.3,  at 
which  the  first  peak  of  the  body  lift  is  observed.  Two  low-pressure  regions  on  the  thorax  and 
posterior  body  can  be  clearly  observed  at  this  moment,  which  coincide  with  the  locations  of 
TXVs  and  PBVs,  implying  that  they  were  generated  by  the  two  vortices.  At  the  time  instant 
t/T=0.74,  at  which  the  second  lift  peak  occurs,  the  low-pressure  regions  generated  by  the  PBVs 
disappears  (Figure  14b),  leaving  only  one  low-pressure  region  on  the  body.  This  is  because 
PBVs  have  already  detached  with  the  body.  In  addition,  during  upstroke,  this  low-pressure 
region  moves  to  the  top  surface  of  the  body  and  follows  the  movement  of  TXVs  (Figure  14c). 
According  to  this  analysis  about  the  spacial  and  temporal  change  of  the  pressure  on  the  body,  the 
high  body  lift  generated  during  downstroke  is  attributed  to  both  TXVs  and  PBVs,  while  the  lift 
peak  which  occurs  at  the  mid-upstroke  is  mainly  due  to  the  TXVs. 


Figure  14.  Surface  pressure  coefficient  distribution  on  the  body  for  the  WB  model  at  selected 
instants:  (a)  t/T=0.3;  (b)  t/T=0.74;  (c)  t/T=0.9. 


Figure  15.  Pressure  difference  between  the  upper  and  lower  wing  surface  (A C  )  of  the  right 
wing  for  the  WB  model  (a)  and  WN  model  (b),  respectively,  (c)  Shows  the  difference  A Cp  of 
between  (a)  and  (b). 
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(a)  (b) 


Figure  16.  Isosurface  of  pressure  coefficient  ( C  =  -0.2 )  at  t/T=0.3  for  (a)  the  WN  model  and  (b) 
the  WB  model.  Wake  schematic  for  (c)  the  WN  model  and  (d)  the  WB  model  at  t/T=0.3. 

To  compare  the  surface  pressure  on  the  wing,  we  plot  the  difference  of  the  pressure 
coefficient  between  the  lower  surface  and  the  upper  surface  (AC  )  in  Figure  15.  Figure  15(a,b) 

indicate  that  the  overall  features  of  the  pressure  on  the  wing  in  WB  and  WN  are  similar.  The 
outer  half  of  the  wing,  especially  on  the  leading-edge  region,  generates  most  of  the  lift  force.  The 
difference  of  A  C  between  WB  and  WN  is  mainly  located  on  the  wing  root  region  (see  Figure 

15c).  An  average  6.5%  enhancement  in  wing  lift  production  can  be  found  in  the  WB  model.  This 
is  attributed  to  the  aforementioned  enhancement  of  the  strength  of  the  LEVs  near  the  wing  root. 

Isosurface  contours  of  a  low-pressure  level  have  been  plotted  in  Figure  16(a)  and  (b)  for  WN 
and  WB,  respectively,  at  t/T=0.3.  Noticeable  in  both  plots  is  a  large  region  of  low  pressure  right 
behind  the  wings  which  is  due  to  the  LEVs,  RVs,  forewing  tip  vortices  (FTVs)  and  hindwing  tip 
vortices  (HTVs).  However,  for  the  WB  model,  an  additional  low-pressure  region  can  be  found 
behind  the  head  and  the  thorax,  which  extends  to  the  region  between  the  sides  of  the  abdomen 
and  the  hindwings.  Figure  16(c)  and  (d)  show  the  schematic  of  the  correlation  between  vortex 
structures  and  regions  of  low  pressure  in  Figure  16(a)  and  (b),  respectively.  The  major  vortical 
structures  of  the  WN  and  WB  models  are  identified.  According  to  Figure  16,  the  extension  of  the 
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low-pressure  region  is  due  to  the  distinct  and  strong  vortex  structures,  TXVs  and  PBVs,  created 
by  WBI.  As  a  result,  lift  production  on  both  the  body  and  the  wings  were  enhanced.  More  details 
about  this  work  can  be  found  in  Ref.  [22], 

Concluding  remarks 

The  effects  of  wing-body  interaction  (WBI)  on  aerodynamic  performance  and  vortex 
dynamics  have  been  numerically  investigated  in  the  forward  flight  of  cicadas.  Flapping  wing 
kinematics  was  reconstructed  based  on  the  output  of  a  high-speed  camera  system.  Following  the 
reconstruction  of  cicada  flight,  three  models,  wing-body  (WB),  body-only  (BD)  and  wings-only 
(WN),  were  then  developed  and  evaluated  using  an  immersed-boundary-method-based 
incompressible  Navier-Stokes  equations  solver.  Results  have  shown  that  due  to  WBIs,  the  WB 
model  had  a  18.7%  increase  in  total  lift  production  compared  with  the  lift  generated  in  both  the 
BD  and  WN  models,  and  about  65%  of  this  enhancement  was  attributed  to  the  body.  This 
resulted  from  a  dramatic  improvement  of  body  lift  production  from  2%  to  1 1.6%  of  the  total  lift 
produced  by  the  wing-body  system.  Further  analysis  of  the  associated  near-field  and  far-field 
vortex  structures  has  shown  that  this  lift  enhancement  was  attributed  to  the  formation  of  two 
distinct  vortices  shed  from  the  thorax  and  the  posterior  of  the  insect,  respectively,  and  their 
interactions  with  the  flapping  wings.  Simulations  are  also  used  to  examine  the  new  lift 
enhancement  mechanism  over  a  range  of  minimum  wing-body  distances,  reduced  frequencies 
and  body  inclination  angles.  This  work  provides  a  new  physical  insight  into  the  understanding  of 
the  body-involved  lift-enhancement  mechanism  in  insect  forward  flight. 

(5)  Aerodynamic  effects  of  morphing  wings  in  dragonfly  forward  flight 

Although  some  encouraging  progress  have  been  made  in  advance  our  knowledge  on  flapping 
morphing  wings  [23-25],  most  of  the  previous  studies  are  more  focused  on  the  aerodynamic 
performance  rather  than  the  vortical  structures.  So  far,  there  is  still  a  lack  of  qualitative  and 
quantitative  descriptions  on  the  vortex  formation  of  deformable  wings,  which  can  lead  to 
improved  models  for  the  design  of  biomimetic  propulsors,  and  also  provide  a  better  understand 
of  vorticity  transport  mechanisms  of  morphing  wings  in  nature.  The  present  work  is  meant  to  fill 
some  of  the  knowledge  gaps  in  this  regard.  The  high-speed  photogrammetry  system,  3-D  surface 
reconstruction  technology  and  numerical  simulations  are  used  to  reveal  the  effects  of  morphing 
wings  of  a  forward  flight  dragonfly.  Specifically,  the  flapping  morphing  wing  kinematics  of  a 
free-flight  dragonfly  are  measured  and  quantified  first  (see  Figure  17).  We  then  use  the 
reconstructed  model  (original  case)  to  explore  the  effects  of  morphing  wings,  first  by  removing 
camber  while  keeping  the  same  time-varying  twist  distribution  (twist-only  case),  and  second  by 
removing  both  the  camber  and  the  spanwise  twist  (rigid  case).  Numerical  simulations  are  carried 
out  using  an  in-house  immersed-boundary-method-based  direct  numerical  simulation  solver.  To 
get  a  better  understand  of  the  aerodynamic  roles  of  morphing  wings,  the  leading-edge  vortex,  the 
wing  surface  pressure  distribution,  and  wake  structures  were  analyzed  and  compared  in  detail  for 
the  model  wings. 
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(a)  y  (b) 


Figure  17.  (a)  Chord-lines  of  dragonfly  forewing  (dashed  line)  and  hindwing  (solid  line)  at  a  few 
instances,  (b)  Definition  of  wing  orientation  angles,  wing  stroke  angle  (t//)  ,  wing  deviation  angle 
((/)) ,  wing  pitch  angle  (0) ,  and  camber  deformation  (h/  c) .  (c-e)  Time  course  of  wing  kinematics 
of  dragonfly  hindwing,  (c)  Wing  stroke  angle  (i//)  and  wing  deviation  angle  ((f)) .  (d)  and  (e)  are 
the  wing  pitch  angles  ( 9 )  and  camber  deromation  (h/  c)  of  different  sections  along  the  wing  span, 
respectively.  The  downstroke  period  is  shaded  as  a  gray  color. 

Using  the  above-reconstructed  wing  kinematics,  the  aerodynamic  forces  and  power  required 
for  the  flapping  motion  of  three  hindwing  models  were  computed  using  the  CFD  solver,  as 
shown  in  Figure  18.  The  stroke  plane  angle  of  hindwings  were  77°  with  respect  to  the  horizontal 
plane,  which  indicates  that  the  hindwings  flapped  with  nearly  a  vertical  stroke.  Therefore,  the 
positive  vertical  force  was  typically  created  during  the  downstroke,  and  slight  negative  force  was 
formed  during  the  upstroke.  The  magnitude  of  the  instantaneous  vertical  of  the  twist-only  wing 
was  smaller  than  that  of  the  rigid  wing,  as  clearly  shown  in  the  middle  position  of  the  stroke 
when  the  twist  angle  of  the  hindwing  was  largest.  The  instantaneous  vertical  force  of  the  original 
wing  with  camber  variation  was  similar  to  that  of  the  twist-only  wing  during  the  downstroke. 
Due  to  the  cruising  motion,  the  force  component  in  the  horizontal  direction  is  close  to  zero,  and 
only  slight  thrust  was  produced  for  compensating  the  body  fraction  drag  during  the  upstroke.  The 
aerodynamic  power  consumption  in  the  plot  shows  that  the  rigid  wing  requires  more  energy  for 
flapping  in  comparison  with  the  other  models. 
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Figure  18.  Time  course  of  lift  (a),  thrust  (b)  and  power  (c)  of  the  hindwing. 


Concluding  remark 

More  details  about  this  work  can  be  found  in  Ref.  [26].  Major  contributions  are  summarized 
as  follows.  In  this  study,  a  high-speed  photogrammetry  system,  3-D  surface  reconstruction 
technology  and  numerical  simulations  are  used  to  reveal  the  effects  of  morphing  wings  of  a 
cruising  dragonfly.  Specifically,  the  flapping  morphing  wing  kinematics  of  a  free-flight 
dragonfly  are  measured  and  quantified  first.  We  then  used  the  reconstructed  model  to  explore  the 
effects  of  morphing  wings,  first  by  removing  camber  while  keeping  the  same  time-varying  twist 
distribution,  and  second  by  removing  both  the  camber  and  the  spanwise  twist.  Our  simulation 
results  revealed  that  the  surface  deformation  can  improve  the  aerodynamic  functions  in  two  ways: 
1)  improving  the  power  economy  by  preventing  the  tip  vortex  bursting;  and  2)  improving  the 
leading-edge  vortex  attachment  by  suppressing  the  generation  of  the  secondary  vortex.  As  a 
result,  the  spanwise  twist  can  boost  the  aerodynamic  efficiency  up  to  20%,  especially  during  the 
wing  translational  phase. 


(6)  Aerodynamics  of  a  dragonfly  turning  maneuver 

Clever  maneuvers  can  be  commonly  observed  in  insect  flight  for  capturing  food  and/or 
avoiding  predators.  Unlike  most  other  insects  such  as  flies,  wasps  and  cicadas  either  reduced 
their  hind-wings  or  mechanically  coupled  fore  and  hind  wings,  dragonflies  have  maintained  two 
independent-controlled  pairs  of  wings  throughout  their  evolution  [27].  Their  neuromuscular 
system  allows  them  to  actively  change  many  aspects  of  wing  motion  in  a  single  wing,  such  as  the 
angle  of  attack,  stroke  amplitude,  and  stroke  plane.  Most  previous  studies  of  are  focused  on  the 
aerodynamics  of  dragonfly-like  tandem  wings  in  steady  flight  motion.  Although  the  unsteady 
free  flights  of  dragonflies  have  also  been  studies  [28-30],  their  works  were  limited  on  reporting 
the  wing  kinematics  and  associated  flight  dynamics.  The  present  effort  is  meant  to  fill  some  of 
the  knowledge  gaps  in  this  regard.  Specially,  a  high-speed  photogrammetry  system  and  3-D 
surface  reconstruction  technology  [1]  are  used  to  reveal  dragonfly  wing  kinematics  during  a 
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turning  maneuver  flight  (see  Figure  19).  The  aerodynamic  performance  and  the  flow  structures 
(see  Figure  20)  is  then  studied  using  an  in-house  immersed-boundary-method-based 
computational  fluid  dynamics  (CFD)  solver.  This  work  aims  to  investigate  the  aerodynamic 
characteristics  of  forewings  and  hindwings  and  its  associated  forewing-hindwing  interaction 
effects  in  a  turning  maneuvering  motion. 


Figure  19.  Motion  reconstruction  of  dragonfly  taking-off  maneuver.  The  side  panels  show  4  of 
110  frames  recorded  by  high-speed  videography. 

More  details  about  this  work  can  be  found  in  Ref.  [26].  Major  contributions  are  summarized 
as  follows.  The  asymmetrical  wing  kinematics  generate  unbalanced  forces  in  both  vertical  and 
horizontal  directions.  During  the  downstroke,  the  force  generated  by  the  outer  wings  are  higher 
than  inner  wings.  By  contrast,  the  inner  wings  generate  higher  forces  during  the  upstroke.  The 
surface  pressure  distribution  shows  that  the  majority  of  the  force  is  generated  around  the  leading- 
edge.  The  unbalanced  forces  lead  to  the  torques  along  its  body  axis  in  three  directions. 
Specifically,  both  forewings  and  hindwing  contribute  to  the  roll  torque  generation,  whereas 
forewings  are  dominant  in  generating  pitch  and  yaw  torques.  To  examine  the  effect  of  forewing¬ 
hindwing  interaction  on  the  aerodynamic  performance,  two  additional  cases,  one  with  forewings 
only  and  the  other  with  hindwings  only,  are  carried  out.  Our  results  indicate  that  the  interaction 
can  lead  to  a  slight  forces  reduction,  but  can  save  the  aerodynamic  power  up  to  11%. 
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Figure  20.  3-D  vortex  structures  in  the  flow  during  the  early  turning  motion,  where  the  time 
stamp  from  (a)  to  (d)  is  23,  29,  35  and  40  ms.  The  vortex  loop  from  the  downstroke  is  marked  by 
dashed. 


2.3  Low-dimensional  Modeling  and  Analysis  of  the  Morphing  Wing  Kinematics 

1}  Spherical-coordinates-based  Singular  Value  Decomposition  [SSVDj 

Singular  value  decomposition  (SVD,  also  known  as  POD  or  PCA  in  some  fields  of 
application)  is  a  powerful  method  for  data  analysis  aimed  at  obtaining  low-dimensional 
approximate  descriptions  of  a  high-dimensional  process  or  dataset  [31].  The  most  remarkable 
feature  of  the  SVD  is  its  optimality:  it  provides  the  most  efficient  way  of  capturing  the  dominant 
components  of  any  dataset  with  only  a  finite  and  often  surprisingly  few  numbers  of  modes.  In 
gait  analysis,  PCA  has  yielded  insights  into  human  walking  strategies  and  the  interrelationships 
in  terms  of  temporal,  kinematic  and  kinetic  variables.  Urtasun  and  colleagues  [32]  have  used 
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PCA  to  identify  invariant  or  common  features  within  the  whole  body  kinematics  of  a 
contemporary  dance  movement  pattern.  For  animals  in  nature,  Bozkurttas  and  colleagues  [33] 
have  used  SVD  to  study  the  pectoral  fin  kinematics  and  its  associated  aerodynamics  of  bluegill 
sunfish.  Representing  the  motions  as  a  linear  sum  of  principal  components  has  become  a  widely 
accepted  animation  technique  [34,  35].  SVD  and  other  similar  methods  are  closely  related,  and 
the  close  connections  and  equivalence  of  these  various  methods  can  be  found  elsewhere  [31]. 

SVD  can  be  considered  as  an  extension  of  the  traditional  eigenvalue  decomposition  for  the 
non-square  matrix,  which  contains  dataset  that  represent  the  wing  motion  in  both  time  and  space. 
Displacements  of  all  m  nodes  on  the  wing  surface  at  n  distinct  instants  in  time  are  stored  in  this 
matrix,  named  displacement  matrix.  The  displacement  matrix  (denoted  by  A )  is  as  follows: 
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The  displacements  stored  in  above  matrix  are  calculated  as  follows: 

A  r‘=r‘-r^ 

<A  dt=G,-OreS  (9) 

a 

where  denote  the  coordinates  of  the  node  5  at  time  instant  t .  Note  that  the  spherical 

coordinate  system  is  used  here.  [r"ef  ,6"ef  ,(f>"ef )  denote  the  coordinates  of  the  node  i  at  a  specific 


reference  time  instant.  An  SVD  of  the  displacement  matrix  A  can  then  be  factorized  as 
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where  U3mx3m  and  Vj/H  are  two  orthogonal  unitary  matrices;  S3mxn  is  a  diagonal  matrix  in  which 
the  diagonal  values  are  called  the  singular  values  of  A ,  which  are  unique.  The  diagonal  elements 
S;i  consist  of  r  =  min(m,  n)  non-negative  numbers  At  ,  which  are  arranged  in  descending  order, 


i.e.  \>A1>--->An>0.  Within  the  SSVD  procedure,  the  A/  values  are  the  square  roots  of  the 

eigenvalues  of  AAr ,  whereas  the  eigenvectors  of  AAr  make  up  the  columns  of  U  and  Vr 
respectively.  In  the  above  expression,  V  represents  the  change  in  each  mode  with  time,  and  U 
contains  the  eigenvectors  corresponding  to  the  spatial  distribution  of  the  modes.  The  singular 
values  A.  can  be  interpreted  as  the  weight  contributions  of  each  mode  in  the  SSVD.  Thus,  the 
‘shape’  of  any  particular  mode  (say  the  kth  mode)  can  be  extracted  by  zeroing  out  all  the  singular 
values  except  for  the  kth  value,  and  reconstructing  from  the  SVD  as  in  Eqn(10).  Similarly,  lower 
dimensional  (say  rank  K<n  )  approximations  to  the  dataset  can  be  obtained  by  using  an 
approximation  to  S  denoted  by  SA.  wherein  AK+l  =  AK+2  =  ■  •  •  =  An  =  0  and  reconstructing  from  the 
SSVD  as  follows: 


A,=US,Vr 


(11) 


2 )  Low -dimensional  Morphology  Analysis  of  Flapping  Wing  Flight  in  Nature 
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In  this  section,  we  will  first  show  details  of  the  SSVD  for  the  fore  wing  motion  of  a  hovering 
dragonfly.  The  displacement  matrix  A  is  now  subjected  to  SSVD.  In  the  following  part,  the 
forewing  motion  of  a  hovering  dragonfly  will  be  used  to  demonstrate  the  decomposition  results. 
30  snapshots  of  the  wing  motion  for  one  complete  flapping  cycle  are  used  in  the  decomposition. 
As  expected,  the  SSVD  leads  to  30  distinct  singular  values,  and  the  spectrum  for  the  first  ten 
singular  values  of  the  wing  motion  is  shown  in  Figure  21  along  with  a  cumulative  plot  for  the 

same  data.  The  normalized  singular  value  for  kth  mode  A,.  is  defined  as: 

K  =  (12) 

i= 1 

The  singular  values  are  normalized  by  the  sum  of  all  singular  values.  Therefore,  the 
cumulative  values  sum  to  unity.  A  number  of  interesting  observations  can  be  made  from  this  plot. 
First,  the  singular  value  spectrum  shows  three  distinct  ranges:  the  first  between  mode  1  -  2,  in 
which  we  see  a  rapid  decrease  in  the  amplitude,  the  second  from  mode  2  -  4  in  which  there  is  a 
much  slower  reduction  in  amplitude  and,  finally,  the  range  from  mode  4  -  30  that  has  negligible 
(less  than  2%)  total  contribution.  The  rapid  initial  decrease  in  the  spectrum  is  significant  which 
suggests  that  a  small  number  of  modes  contain  most  of  the  essential  features  of  the  wing  gait.  In 
fact,  the  cumulative  values  show  that  the  first  two  and  three  modes  capture  about  93.1%  and  96.0% 
respectively  of  the  total  motion.  In  fact,  only  the  first  mode  captures  close  to  84.7%  of  the 
motion  of  the  wing,  which  is  a  clear  demonstration  of  the  ability  of  SSVD  to  represent  the 
dataset  with  the  least  possible  number  of  modes. 

The  gait  corresponding  to  individual  modes  can  be  extracted  as  described  above,  and  the 
surface  conformations  for  each  of  these  extracted  modes  are  then  constructed  using  the  original 
wing  mesh  with  triangular  elements.  The  first  two  modes  are  highly  distinct  and  relatively  easy 
to  interpret,  and  we  briefly  describe  the  key  qualitative  features  of  these  modes.  Figure  22(b,  c) 
shows  mode  1  and  mode  2  at  five  different  times  during  one  flapping  cycle.  Also  shown  on  the 
left  for  direct  comparison  are  the  wing  motion  from  the  experiment  (also  called  the  ‘mode-all’ 
case,  since  it  contains  all  the  SSVD  modes).  In  these  figures,  the  colors  reflect  wing  deformation 
by  plotting  contours  of  distances  between  vertices  on  the  wing  surface  to  the  corresponding  least 
square  plane  of  the  wing.  Mode  1  involves  very  large  rotating  motion  about  the  wing  root,  which 
is  called  the  ‘flapping’  motion.  The  wing  flaps  back  and  forth  with  certain  offset  angle.  And  also, 
the  mode  one  shows  minimum  deformation  during  the  stroke.  This  mode  is  actively  produced  by 
the  dragonfly  through  flight  mussels  at  the  wing  root. 
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Figure  21.  The  SSVD  spectrum  of  the  first  ten  modes  for  the  dragonfly  right  forewing  kinematics. 
The  left  ordinate  shows  |/t| / ^|/L|  >  and  the  right  ordinate  shows  the  cumulative  value  of  the  left 
ordinate. 


Figure  22.  Wing  motions  of  (a)  experimental  kinematics  (also  called  the  mode-all  case),  (b) 
SSVD  mode  1,  (c)  mode  2  and  (d)  low-dimensional  model  mode  1+2.  The  wings  are  colored 
with  distances  between  wing  surfaces  and  corresponding  least  square  planes.  The  distances  are 
normalized  by  wing  mid  chord  length. 

Mode  2  is  a  twisted  motion  primarily  in  the  span- wise  direction,  which  occurs  along  the  span 
axis  of  the  wing.  It  presents  as  chord-wise  rotations  of  the  wing  during  the  reversal  phase.  In 
contrast  to  mode  1,  this  mode  is  primarily  a  result  of  flow-induced  deformation.  It  can  be 
deduced  from  the  fact  that  there  are  no  muscles  in  the  wing  surface  that  could  produce 
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deformation  in  the  wing.  Furthermore,  the  deformation  is  primarily  in  the  direction  of  the  flow 
relative  to  the  wing  motion,  which  supports  the  assertion  that  this  mode  is  flow-induced.  The  rest 
of  the  modes  in  the  spectrum  are  associated  with  relatively  small  motions  that  are  not  very 
distinct.  We,  therefore,  do  not  describe  these  individually,  although  we  will  consider  the  effect  of 
mode  3  on  the  aerodynamics  in  the  following  sections. 

SSVD  has  decomposed  the  wing  motion  into  its  orthogonal  components  and  helped  us 
understand  the  key  features  of  the  nature  flyers  wing  motion.  The  SSVD  results  can  also  be  used 
to  reconstruct  low-dimensional  approximations  of  the  mode-all  case  using  a  subset  of  the 
orthogonal  modes.  Lower  dimensional  models  of  the  wing  gait  are  synthesized  by  successively 
adding  modes  to  Mode  1 .  The  forewing  motion  of  a  hovering  dragonfly  is  used  to  demonstrate 
the  method.  Figure  22(a,  b,  d)  shows  the  surface  snapshots  at  four  different  time  instances  during 
one  flapping  cycle  for  mode  1  and  mode  1+2  in  comparison  to  the  complete  (mode-all)  motion. 
The  similarity  between  the  wing  shapes  for  Mode  1+2  and  the  Mode- All/experiment  cases  is 
evident  in  this  figure.  Removal  of  higher  SSVD  modes  from  the  motion  is  analogous  to  filtering 
the  experimental  data  in  space  and  time. 

The  SSVD  analysis  suggests  one  natural  approach  to  the  development  of  the  robotic  wing. 
Since  a  small  number  of  modes  capture  a  significant  portion  of  the  motion,  it  stands  to  reason 
that  a  systematic  procedure  for  developing  a  robotic  wing  would  involve  designing  actuation 
mechanisms  that  reproduce  a  small  number  of  these  modes.  The  question  that  remains  to  be 
answered  is  what  kind  of  performance  can  we  expect  from  these  lower  dimensional  wing  models, 
and  how  does  the  performance  scale  as  we  include  additional  modes?  This  will  allow  us  to  make 
a  rational  compromise  between  the  complexity  of  wing  design  and  wing  performance.  It  should 
be  noted  here  that  the  performance  is  a  consequence  of  the  flow  associated  with  these  low¬ 
dimensional  wing  models.  Thus,  even  though  the  modes  are  kinematically  linear  (and  therefore 
additive),  the  performances  are  not  expected  to  scale  linearly  with  the  modes  since  the  flow  is 
governed  by  the  Navier-Stokes  equations  that  are  nonlinear.  Thus,  the  answer  to  the  above 
question  requires  that  we  explicitly  determine  the  performances  for  these  low -dimensional  wing 
models.  The  following  sections  describes  our  approach  to  answering  this  question. 

We  further  study  the  effects  of  increasing  the  dimensionality  of  the  fore  wing  motion  of  the 
hovering  dragonfly  on  the  aerodynamic  performances.  The  effects  of  model  dimensionality  on 
the  quantitative  characteristics  of  the  wing  are  investigated,  including  force  production  and  lift 
economy. 

The  time  variations  of  lift,  drag  and  power  coefficients  are  presented  for  all  the  low¬ 
dimensional  gaits  and  compared  to  the  mode-all  case  in  Figure  23(a  -  c),  respectively.  Several 
observations  on  how  each  SSVD  mode  contributes  to  the  performances  of  the  wing  can  be  made 
from  these  results.  It  should  be  noted  that  only  mode  1  can  be  simulated  by  itself.  However, 
given  the  underlying  nonlinearity  of  the  flow,  the  contribution  of  mode  2  and  mode  3  are 
investigated  by  considering  the  differences  in  the  performances  from  the  lower-level  gait.  Thus, 
the  effect  of  mode  2  on  performance  is  obtained  by  analyzing  the  differences  between  the 
performances  of  the  mode  1  and  mode  1+2  cases.  Similarly,  the  effect  of  mode  3  on  wing 
performances  can  be  assessed  by  comparing  the  performances  of  the  mode  1+2+3  case  with  that 
of  the  mode  1+2  case. 
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(a)  (b)  (c) 


Figure  23.  Comparison  of  the  time  variation  of  aerodynamic  performances  between  the  mode-all 
and  low-dimensional  gaits,  (a)  lift;  (b)  drag;  (c)  power. 

Table  2.  Cycle  averaged  aerodynamic  performance  of  the  mode-all  and  low-dimensional 

gaits. 


r(%) 

Q 

c 

^ PW 

7 

Lift  production 

(%) 

Mode-all 

100 

0.752 

1.021 

0.776 

0.969 

100 

Mode  1 

84.7 

0.149 

1.464 

1.087 

0.137 

19.8 

Mode  1+2 

93.1 

0.749 

1.100 

0.818 

0.916 

99.6 

Mode  1+2+3 

96.0 

0.750 

1.082 

0.813 

0.923 

99.7 

As  we  can  see  in  Figure  23  that  the  aerodynamic  performances  are  very  similar  for  all  cases 
except  for  the  mode  1  only  case.  For  the  lift  production.  Figure  23(a)  shows  that  all  cases 
produce  positive  lift  except  for  mode  1  case,  which  negative  lift  can  be  observed  during  the 
upstroke.  Also,  the  lift  production  for  mode  1  case  during  downstroke  is  much  smaller  (about  2.5 
times  smaller)  than  other  cases.  For  the  drag  production,  Figure  23(b)  shows  that  the  drag  is 
much  higher  for  mode  1  only  case  during  the  downstroke.  For  upstroke,  the  drag  productions  are 
very  similar  for  all  cases.  These  are  because  the  lack  of  wing  rotation  about  the  span  axis,  which 
is  included  in  the  deformation  mode  (mode  2).  During  the  downstroke,  the  wing  angle  of  attack 
for  mode  1  case  is  much  greater  than  other  cases,  and  it  leads  to  less  lift  and  more  drag 
production  (Figure  23b).  During  the  upstroke,  the  wing  angle  of  attack  for  mode  1  case  is  greater 
than  90  degrees  due  to  lack  of  wing  rotation  about  the  span  axis,  and  negative  lift  produces  at 
this  phase.  The  involvement  of  the  deformation  mode  (mode  2)  can  greatly  improve  the  lift 
production  and  reduce  the  drag  produced  by  the  wing. 

For  the  aerodynamic  power  histories  (Figure  23c),  we  can  see  that  all  cases  present  two  peaks 
during  the  cycle,  and  the  amplitude  of  power  consumption  for  mode  1  case  is  much  greater  than 
other  cases  at  both  down  and  upstrokes.  This  suggests  that  involving  the  deformation  mode 
(mode  2)  can  also  reduce  the  power  consumption  of  the  dragonfly. 

Cycle  averaged  aerodynamic  performances  are  listed  in  Table  2.  We  can  see  from  the  table 
that  mode  1+2  is  a  good  approximation  of  the  original  wing  motions  mode-all.  It  contains  only 
two  dominant  SSVD  modes,  and  the  motion  is  recovered  over  90%.  The  associate  aerodynamic 
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performances  of  mode  1+2  case  are  very  similar  to  the  mode-all  case  as.  The  lift  production  is 
recovered  93%  of  the  original  motion. 

For  the  wake  structures,  we  cut  slices  along  the  wingspan  to  see  the  leading  edge  vortex 
structures  and  measure  the  associated  circulations  evaluated  at  each  slice.  Figure  24  shows  the 
corresponding  results.  We  can  see  from  the  figure  that  for  all  cases,  the  LEV  shapes  gradually 
grow  bigger  from  the  wing  root  to  wing  tip.  More  importantly,  the  corresponding  vortex 
structures  are  very  different  in  mode  1  only  case  comparing  to  other  cases.  The  LEV  shapes  are 
much  bigger,  and  the  associated  attachment  is  bad  in  mode  1  case.  For  other  cases,  the  LEV 
shapes  are  similar.  Small  differences  can  be  observed  near  the  wing  tip  region.  At  mid 
downstroke  (t/T=0.25),  the  LEV  attachment  is  pretty  good  for  all  cases  except  for  the  mode  1 
only  case.  However,  at  mid  upstroke  (t/T=0.75),  the  LEV  attachment  is  not  as  good  as  that  at  mid 
downstroke. 


t/T= 0 


75 


mode  1+2 


t/T=  0.25 


(b) 


mode-all 


(d) 


mode  1 


Figure  24.  2D  slices  of  the  leading  edge  vortex  along  the  wingspan  (10%~90%)  at  two  time 
instances  for  (a,  b)  mode-all,  (c,  d)  mode  1,  and  (e,  f)  mode  1+2  cases.  The  corresponding  vortex 
center  are  marked  with  green  dots  at  each  slice.  The  contours  represent  normalized  span-wise 
vorticity. 
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Quantitative  measurements  related  to  the  LEV  attachments  of  all  cases  discussed  above  are 
conducted.  We  first  determine  the  LEV  centers  based  on  vortex  shapes  shown  in  Figure  24.  After 
that,  we  measure  the  distances,  which  are  named  as  lift-off  distances,  between  LEV  centers 
(green  dots)  and  the  wing  surface  to  evaluate  the  LEV  attachments.  Figure  25  shows  the  results 
at  mid  downstroke  (t/T=0.25)  and  mid  upstroke  (t/T=0.75).  For  the  mode  1  only  case,  the  lift-off 
distances  are  much  higher  than  other  cases  in  both  time  instances,  which  indicates  bad  LEV 
attachment.  In  addition,  the  ranges  of  the  lift-off  distance  are  much  wider  in  mode  1  only  case  for 
both  time  instances.  It  ranges  from  0.15  chords  to  0.55  chords  at  mid  downstroke,  and  from  0.20 
chords  to  1.09  chords  at  mid  upstroke.  Moreover,  at  mid  downstroke,  a  small  peak  of  lift-off 
distance  shows  up  at  60%  span  and  then  drops  at  70%  span,  which  corresponds  to  the  shed  of 
LEV  at  70%  span  in  Figure  24(c). 

For  the  mode-all  and  mode  1+2  cases,  the  lift-off  distances  are  quite  similar.  The  differences 
are  less  than  8%  and  10%  for  each  span  location  at  mid  downstroke  and  mid  upstroke, 
respectively.  More  importantly,  at  both  time  instances,  two  distinct  ranges  of  the  lift-off 
distances  can  be  observed.  The  first  range  is  from  10%  to  60%  span,  where  the  lift-off  distance 
increases  slowly.  It  suggests  that  the  LEV  attached  pretty  well  with  this  range.  The  second  range 
is  from  60%  to  90%  span,  where  a  rapid  increase  in  lift-off  distance  can  be  observed.  The  LEV  is 
lifted  by  the  tip  vortex  and  starts  to  merges  with  the  tip  vortex  at  this  point. 


Span  (%)  Span  (%) 

Figure  25.  Distances  of  LEV  center  to  the  wing  surface  (lift-off  distances)  at  (a)  t/T=0.25  and  (b) 
t/T=0.75  for  mode-all,  mode  1,  and  mode  1+2  cases. 


The  lift-off  distances  at  mid  downstroke  are  almost  twice  as  much  as  that  at  mid  upstroke  for 
all  cases  and  span  locations,  which  indicates  that  the  LEV  attachment  at  downstroke  is  much 
better  than  upstroke. 

Quantitative  measurements  of  LEV  circulation  distributions  along  the  wingspan  are  also 
performed  for  all  cases  discussed  above  based  on  the  2D  flow  slices  shown  in  Figure  24.  Figure 
26  shows  the  corresponding  results  at  mid  downstroke  (t/T=0.25)  and  mid  upstroke  (t/T=0.75). 
The  circulation  is  calculated  and  normalized  as  follows: 


r*  = 


(13) 


Where  S  stands  the  surface  of  LEV  shapes  shown  in  Figure  24;  ®  is  vorticity  on  S  ;  Uref  is 


reference  velocity,  which  is  chosen  as  the  average  velocity  of  wing  mid  chord;  c  denotes  the  mid 
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chord  length.  We  can  see  from  Figure  26  that  the  LEV  circulations  of  mode-all,  mode  1+2  and 
mode  1+2+3  cases  are  very  close.  The  difference  is  less  than  7%  at  mid  downstroke  and  3%  at 
mid  upstroke.  For  all  cases,  the  LEV  circulations  gradually  increase  from  the  wing  root  to  wing 
tip  and  drops  a  little  bit  near  the  wing  tip  region.  Maximum  circulation  can  be  observed  at  around 
80%  span,  and  the  corresponding  value  for  mode  1  only  case  is  about  1.5  and  1.3  times  larger 
than  that  of  other  three  cases  at  mid  downstroke  and  mid  upstroke,  respectively. 
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Figure  26.  LEV  circulations  at  (a)  t/T=0.25  and  (b)  t/T=0.75  for  mode-all,  mode  1,  and  mode 
1+2  cases. 

In  summary,  we  have  studied  three-dimensional  flows  around  the  forewing  of  a  hovering 
dragonfly  in  this  section,  focusing  on  the  effects  of  model  dimensionality.  Both  qualitative 
observations  and  quantitative  measurements  are  performed.  The  results  show  that  for  cases  of 
mode-all,  mode  1+2  and  mode  1+2+3,  the  wake  structures  in  both  far  and  near  field  are  similar, 
while  significant  differences  can  be  found  in  mode  1  only  case.  Quantitative  measurements  of 
the  flow  field  at  two  time  instances  (t/T=0.25  and  0.75)  are  performed  in  all  cases,  including  the 
LEV  lift-off  distances  and  circulations.  The  results  show  the  similarity  in  all  cases  except  for  the 
mode  1  only  case,  which  has  much  greater  LEV  lift-off  distances  and  circulations. 

Concluding  Remarks: 

In  this  work,  a  combined  experimental  and  computational  method  is  developed  to  study  the 
complex  morphing  wing  kinematics  and  its  associated  aerodynamics  of  a  hovering  dragonfly. 
SSVD  analysis  shows  that  the  wing  motion  can  be  described  by  only  three  dominant  modes.  The 
first  two  dominant  modes,  which  contain  93.1%  of  the  total  motion,  are  highly  distinguishable. 
The  mode  1  (flapping  mode)  consists  of  a  simple  flapping  motion,  and  the  mode  2  (morphing 
mode)  contains  morphing  in  span-wise  and  chord-wise  directions.  The  simulation  results  show 
that  with  the  help  of  the  mode  2,  the  lift  production  and  lift  economy  are  greatly  improved 
comparing  to  the  mode  1  only  case.  Also,  the  mode  1+2  case  can  recover  96%  of  the  lift 
production  and  91%  of  the  lift  economy  of  the  original  case.  By  studying  the  unsteady  flow  field 
of  all  cases,  we  conclude  that  the  SSVD  mode  2  can  greatly  reduce  the  large  tip  vortex  found  in 
the  mode  1  only  case.  In  addition,  the  leading  edge  vortex  attachment  is  greatly  improved  for  the 
cases  with  the  mode  2. 
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2.4  Optimization  of  Low-dimensional  Morphing  Wing  Models 

In  this  section,  we  investigate  the  optimal  configurations  of  dominant  modes  on  aerodynamic 
performance  for  the  dragonfly  wing.  Figure  27  shows  the  convergence  histories  of  the  three 

design  variables,  W1 ,  W2 ,  and  cp  .  The  objective  functions  are  chosen  as  the  cycle  averaged  lift 
coefficient  and  the  lift  efficiency.  The  initial  guesses  of  the  design  variables  correspond  to  the 

low  dimensional  model  Mode  1+2  (W,  =1,  W2  =  1 ,  and  ^  =  0  ).  We  can  see  from  the  figure  that 
the  two  optimization  case  converge  within  6  iterations. 


w,  w2 


) 


Figure  27.  Convergence  history  for  design  variables  (a)  Wj ,  (b)  W2 ,  and  (c)  cp  of  the  two 
optimization  cases. 

Figure  27  shows  the  time  histories  of  aerodynamic  performance  of  the  two  optimization 
cases  along  with  the  case  of  low  dimensional  model  Mode  1+2  for  comparison.  We  can  see  that 
the  case  Opt  CL  shows  largest  amplitude  of  the  lift  production  during  the  downstroke  among  the 
three  cases,  while  the  case  Opt  p  presents  lowest  of  that.  During  the  upstroke,  the  case  Opt  p 
shows  largest  amplitude  of  the  lift  production  among  the  three  cases,  while  the  case  Mode  1+2 
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presents  lowest  of  that.  For  the  power  consumption,  the  case  Opt  rj  shows  the  lowest  power 
consumption. 


t/T  t/T 

Figure  28.  Comparison  of  the  time  course  of  aerodynamic  performance  of  Mode  1+2  and 
optimized  wing  gaits,  (a)  Lift;  (b)  power. 

Table  3  lists  the  cycle  averaged  aerodynamic  performance  for  the  two  optimum  cases  and  the 
case  Mode  1+2.  Comparing  to  the  case  Mode  1+2,  the  thrust  production  of  the  Opt  CL  case  is 
increased  by  7.7%,  and  the  propulsive  efficiency  of  the  Opt  i)  case  is  increased  by  51.6%.  The 
SSVD  modes  can  greatly  improve  the  aerodynamic  performance  of  the  flapping  wing,  especially 
for  the  lift  efficiency. 

Table  3.  Cycle  averaged  hydrodynamic  performance  of  Mode  1+2  and  optimized  wing  gaits. 


Cases 

Q 

r 

^ PW 

i] 

Mode  1+2 

0.750 

0.818 

0.917 

Opt  CL 

0.808 

0.730 

1.107 

Opt  rj 

0.738 

0.531 

1.390 

Figure  29  shows  the  wake  structure  at  the  t/T=0.27  (near  the  mid  downstroke)  and  t/T=0.73 
(near  the  mid  upstroke)  of  the  three  cases.  We  can  see  that  the  general  wake  structures  are  similar 
for  all  the  cases.  However,  the  case  Opt  CL  shows  the  strongest  LEV  and  case  Opt  //  shows  the 
weakest  wing  tip  vortex. 
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Figure  29.  Wake  structures  for  the  (a)  Mode  1+2,  (c)  Opt  CL,  (e)  Opt  rj  at  t/T=0.27;  wake 
structures  for  the  (b)  Mode  1+2,  (d)  Opt  CL ,  (f)  Opt  //  at  t/T=0.73. 
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Concluding  Remarks: 

In  this  work,  we  have  investigated  the  optimal  configurations  of  dominant  SSVD  modes  on 
aerodynamic  performance  for  the  dragonfly  wing.  The  optimized  low  dimensional  wing  models, 
which  can  beyond  biological  levels  of  aerodynamic  performance,  are  obtained.  The  associated 
flow  mechanisms  are  found  to  be  the  improved  LEV  strength  and  the  reduced  TV  strength. 


3  Adjoint  Optimization 


The  mechanism  of  flapping-wing  aerodynamics  provides  an  efficient  way  to  generate  necessary 
lift  and  thrust  and  is  the  most  common  way  of  flying  adopted  by  birds  and  insects.  In  comparison 
to  other  flying  mechanisms  used  by  nature’s  flyers  and  artificial  flying  machines,  flapping  shows 
many  attractive  characteristics  such  as  agility,  hovering  capability,  efficiency  at  low  Reynolds 
number,  etc.  Although  some  understanding  is  achievable  through  carefully  designed  numerical 
simulations  [[36]  Dong,  H.,  Mittal,  R.  &  Najjar,  F.  2006  Wake  topology  and  hydrodynamic 
performance  of  low  aspect-ratio  flapping  foils.  J.  Fluid  Mech.  566,  309-343. 

-40],  the  problem’s  large  parametric  space  prevents  further  physical  understanding  and 
optimization  through  a  direct  parametric  study. 

To  achieve  an  understanding  of  all  the  control  parameters,  one  often  chooses  to  reduce  the 
complexity  of  the  physical  model  and/or  the  size  of  parametric  space.  Based  on  a  quasi-steady 
model  with  11  control  parameters,  Berman  &  Wang  [41]  were  able  to  use  a  hybrid  algorithm  of 
the  genetic  method  and  simplex  method  to  minimize  the  power  consumption  of  insect  flights. 
Ghommem  et  al.  [42]  used  an  unsteady  vortex-lattice  method  and  a  deterministic  global 
optimization  algorithm  for  the  optimization  of  flapping  wings  in  forward  flight  with  active 
morphing,  where  only  4-8  parameters  were  considered.  Milano  &  Gharib  [43]  applied  a  genetic 
algorithm  in  an  experimental  setting  to  maximize  the  average  lift  from  a  flapping  flat  plate  by 
limiting  the  number  of  control  parameters  to  only  4.  Trizila  et  al.  [44]  has  used  a  combined 
approach  with  numerical  simulation  and  surrogate  modelling  to  explore  a  three-parameter  design 
space  for  a  three-dimensional  plate  in  hovering  motion.  Building  a  map  of  the  entire  desire  space 
is  useful  for  some  problems  but  may  not  be  necessary  for  others.  On  the  other  hand,  the 
computational  cost  was  still  high  even  with  a  surrogate  model  and  prevents  the  study  from 
including  more  parameters,  and  the  accuracy  was  limited  by  the  surrogate  model.  There  have 
also  been  efforts  to  use  gradient-based  methods  for  optimization.  Tuncer  &  Kaya  [45]  and 
Culbreth,  Allaneau  &  Jameson  [46]  used  a  gradient-based  method  to  optimize  the  flapping- wing 
motion  for  better  thrust  and  efficiency.  They  used  direct  numerical  simulation  for  each  set  of 
control  parameters  and  then  computed  the  gradient  of  the  cost  function  subject  to  the 
perturbation  of  parameters  directly  by  finite  difference.  Since  all  parameters  were  evaluated 
individually  and  required  their  own  simulation,  the  process  became  extremely  expensive,  and  the 
former  work  was  limited  to  4  parameters  and  the  latter  included  cases  of  from  1  to  1 1  parameters. 

Different  from  the  above  works,  an  adjoint-based  method  is  capable  of  obtaining  the  gradient 
information  simultaneously  for  an  arbitrary  number  of  input  parameters  by  one  single 
computation  in  the  adjoint  space.  Consequently,  the  total  computational  cost  to  obtain  the 
sensitivity  of  a  cost  function  to  all  control  parameters  is  independent  of  the  number  of  control 
parameters.  Thus  an  adjoint-based  method  is  suitable  for  the  sensitivity  analysis  and  optimization 
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of  problems  with  a  large  input  space  but  a  small  output  space,  such  as  the  study  of  kinematics 
and  deformation  of  flexible  flapping  wings.  Depending  on  the  order  of  applying  the 
discretization  and  adjoint  formulation,  there  are  two  types  of  adjoint  approaches:  the  continuous 
approach  [47,48]  and  the  discrete  approach  [49,50],  Nadarajah  &  Jameson  [51]  and  Collis  et  al. 
[52]  have  discussed  the  pros  and  cons  of  the  two  approaches.  In  the  current  work,  we  take  the 
continuous  approach  for  its  simplicity  and  clarity  in  the  governing  equation  for  adjoint  space, 
which  has  mathematical  terms  which  may  be  interpreted  as  representing  the  generation, 
convection  and  dissipation  of  adjoint  variables  [53],  However,  use  of  the  continuous  adjoint 
approach  on  flapping  wings  has  been  scarce  because  of  the  difficulty  in  defining  the  perturbation 
at  a  moving  or  deforming  boundary  [54],  In  their  work  of  applying  an  adjoint-based  method  to 
obtain  the  gradient  information  for  the  shape  optimization  of  a  plunging  airfoil,  Nadarajah  & 
Jameson  [55]  used  a  mapping  function  to  transfer  the  physical  domain  with  a  moving  boundary 
to  a  computational  domain  with  a  fixed  boundary  so  that  traditional  adjoint-based  methods  for  a 
fixed  domain  can  be  directly  applied  (and  the  trouble  relating  to  the  moving  boundary  is 
avoided).  Although  the  idea  was  straightforward,  the  mapping  function  increased  considerably 
the  complexity  of  the  formulation  even  for  their  simple  case  where  only  the  shape  (as  a  steady 
function)  of  a  rigid  flapping  airfoil  was  optimized.  In  other  words,  Nadarajah  and  Jameson  were 
only  able  to  optimize  the  steady  part  of  an  unsteady  mapping  function.  When  the  moving 
trajectory  and  dynamics  morphing  (as  an  unsteady  function)  need  to  be  optimized  for  a  flapping 
wing,  the  complexity  reaches  a  much  higher  level  and  yields  the  mapping-function  approach 
infeasible.  To  deal  with  the  increased  complexity  introduced  by  an  unsteady  morphing  domain, 
we  choose  to  apply  non-cylindrical  calculus  [54,56]  to  derive  adjoint  equations  directly  in  a 
morphing  domain  and  to  optimize  the  moving  boundary  in  its  original  space  without  using  a 
mapping  function.  The  advantage  of  choosing  non-cylindrical  calculus  over  the  unsteady 
mapping  function  has  been  demonstrated  by  Protas  &  Liao  [56]  with  a  comparison  of  the  two 
approaches  in  deriving  the  adjoint  equation  for  a  one-dimensional  heat  equation  with  a  moving 
boundary. 

Using  a  continuous  adjoint  approach  to  handle  the  large  parametric  space  and  non-cylindrical 
calculus  to  handle  the  moving  boundary,  we  can  study  the  optimal  moving  trajectory  and 
arbitrary  deformation  of  a  flapping  wing.  The  optimal  solutions  of  different  configurations  (a 
two-dimensional  rigid  flapping  plate,  a  two-dimensional  flexible  flapping  plate,  a  three- 
dimensional  rigid  flapping  plate  and  a  three-dimensional  flexible  flapping  plate)  for  different 
control  goals  (thrust  performance,  propulsive  efficiency  and  lift  performance  in  hovering) 
provide  a  unique  opportunity  to  understand  the  flapping-wing  mechanism  and  the  role  of 
flexibility  through  a  detailed  comparison  of  optimal  and  non-optimal  controls,  corresponding 
flow  fields  and  other  aerodynamic  performance  indicators. 

3.1  Governing  equations  for  adjoint  optimization 

The  basic  derivation  and  notation  of  the  continuous  adjoint  equations  are  similar  to  those 
used  by  Bewley  et  al.  [47]  and  Wei  &  Freund  [48],  although  those  works  dealt  only  with 
problems  with  a  fixed  domain.  The  inclusion  of  non-cylindrical  calculus  to  formulate  the  adjoint 
equations  a  with  morphing  domain  or  moving  boundary  has  been  suggested  by  Moubachir  & 
Zolesio  [54]  within  a  mathematical  framework  only  and  later  by  Protas  &  Liao  [56]  with 
numerical  implementation  for  simple  problems  (e.g.  the  one-dimensional  heat  equation).  We 
followed  the  same  idea  and  extended  its  application  to  work  for  the  Navier-Stokes  equations. 
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1 )  Governing  equation  and  cost  function 


Though  the  approach  is  not  limited  to  a  particular  shape  or  motion,  for  the  convenience  of 
discussion,  we  consider  a  scenario  where  a  plate  is  plunging  and/or  pitching  with  a  prescribed 
velocity  V(t).  The  flow  dynamics  is  described  by  the  incompressible  Navier-Stokes  equations, 

N{q)  =  F  inQ, 

u  =  V  on  5,  (14) 

uL=uo  in  a 

where  fl  denotes  the  fluid  domain,  S  denotes  the  solid  boundary,  the  primary  variable  q=[u,  p] 
for  pressure  and  velocity,  the  operator 

dit 
dx, 


N(q)  = 


du;  d(ujUj) 

dt  dx. 


d  u.  dp 
V - T  +  — 

dx,  dx, 


F  =  0, 


j  j 

For  demonstration  purpose,  we  pick  a  simple  cost  function  J,  which  is  to  minimize  the 
overall  difference  between  the  velocity  u  at  a  downstream  region  fl0  and  a  target  velocity  u^oin 
that  region  for  time  period  (0,  T): 


/=J.JJU-U< 


q> 


dCldt 


(16) 


This  choice  of  cost  function  provides  the  convenience  of  having  an  obvious  optimal  solution: 
u  =  ufi((,  which  makes  an  easy  validation  of  the  approach. 

2]_Adjoint  equation  and  gradient  calculation 


Introducing  non-cylindrical  calculus  to  define  the  boundary  perturbation  and  adjoint 
variables  q*=[u*,  p*]  and  Z*are  as  Lagrange  multipliers,  we  get  derive  the  adjoint  equation  with 
controls  for  moving  boundary: 

N*(q)q*  =  F*  inQ, 
u*  =0  on  S  Ur^, 
u  |  =  0  in  Q, 

(17) 


I  t=T 


dZ_ 

dt 


+  Z*divsY  =  (Vu)r  •a*n  on  5, 


I  t=T 


=  0  on  5, 


where 
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N*(q)q  = 


8u*  duj  |  d2u*  dp * 


— -  +  u.  — -  +  — -  +v — ^  +  -!— 

dt  ;  dx ■  dx.  dx2  dx, 


2(ui 


1  in 

0  otherwise’ 


&  =~pS  +u(^  +  —t). 

,J  ,J  dxj  dx, 

The  above  formulation  defines  the  adjoint  equations  with  a  set  of  desirable  boundary  and 
initial  conditions,  where  t=T  is  considered  the  ‘initial’  time,  since  the  adjoint  system  typically 
evolves  backwards  in  time.  With  the  definition  of  the  adjoint  system,  the  gradient  is  given  by: 


J'  =  - 


where 


+z*)dsdt’ 


*  *  /  Xsvlj  1  \ 

an=  p  Sn+u(-^  +  -^). 


When  the  control  is  related  to  the  change  of  boundary  location  and  corresponding  velocity,  it 
provides  a  way  to  compute  gradient  information  as  a  function  of  adjoint  solutions  only. 

It  is  noticed  that  the  first  term  in  the  gradient  is  for  the  perturbation  of  the  boundary  velocity 
and  is  the  same  as  the  term  derived  by  the  traditional  approach  for  a  fixed  domain  (which  can 
still  have  velocity  control  on  the  boundary),  and  the  second  term  is  for  the  domain  variation  and 
is  unique  for  a  moving  boundary/domain  and  the  current  approach. 


3.2  Numerical  algorithm 


We  applied  different  moving-boundary  treatments  to  compute  the  forward  (physical)  flow 
field  and  the  backward  (adjoint)  ‘flow’  field,  including  (i)  immersed  boundary  method  [57],  (ii) 
arbitrary  Lagrangian-Eulerian  method  [58]  and  (iii)  for  simple  rigid  cases,  a  moving  reference 
frame  algorithm  [59].  The  numerical  simulations  using  different  boundary  treatments  yield 
similar  results  and  efficiency. 

A  staggered  Cartesian  mesh  with  local  refinement  through  stretching  functions  is  chosen  for 
the  benefit  of  both  computational  efficiency  and  numerical  stability.  We  apply  a  second-order 
central  difference  scheme  for  spatial  discretization,  a  third-order  Runge-Kutta/Crank-Nicolson 
scheme  for  time  advancement  and  a  typical  projection  method  for  incompressible  flow 
conditions  [60].  With  the  similarity  shown  in  the  form  of  adjoint  equations  and  flow  equations, 
we  apply  similar  numerical  algorithms  to  solve  the  adjoint  equations  backwards  in  time,  which 
leads  to  similar  computational  cost  for  the  adjoint  computation. 

The  gradient  can  be  achieved  after  both  the  flow  and  adjoint  equations  are  solved  once.  In 
some  cases,  when  the  control  space  is  large  and  the  gradient  is  not  smooth,  it  may  be  necessary 
to  apply  a  Sobolev  inner  product  or  Savitzky-Golay  smoothing  filters  to  smooth  out  the  gradient. 
Meanwhile,  in  the  main  iteration  to  update  the  gradient,  the  Polak-Ribiere  variant  of  the 
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conjugate  gradient  method  is  used;  within  each  main  iteration,  Brent’s  method  is  used  to 
determine  the  optimal  step  size  along  each  direction  and  this  process  requires  a  couple  of 
subiterations. 

3.3  Validation  of  adjoint-based  optimization  solution  and  gradient  information 

To  assure  the  accuracy  of  the  current  approach,  in  this  section  we  first  compare  the  optimal 
solution  given  by  the  adjoint-based  algorithm  to  a  known  optimal  solution,  then  make  a  further 
comparison  between  the  gradient  information  computed  by  the  adjoint  approach  and  the  gradient 
computed  by  a  ‘brute-force’  finite  difference  approach. 

Here,  we  consider  a  rigid  plate  plunging  normally  to  an  incoming  flow.  The  control  4>  is  the 
plunging  velocity  at  each  individual  time  moment.  We  pick  the  same  cost  function  as  discussed 
in  section  3.1  for  easy  validation. 


Figure  34.  A  typical  snapshot  showing  the  plunging  plate,  the  flow  and  the  observation  region. 

For  validation,  we  first  perform  a  simulation  while  using  the  target  control  4>0,  and  record  the 
velocity  field  Un0  in  fl0  as  the  target  velocity.  With  some  arbitrary  initial  control  4b ,  the 
velocity  in  is  u.  The  difference  between  u  and  unQ,  as  indicated  by  the  cost  function,  should 
drive  the  control  to  match  4>o  if  the  optimization  algorithm  works. 

Figure  35  shows  the  performance  of  the  optimization  algorithm  by:  (i)  checking  the  change 
of  the  cost  function  and  the  gradient  norm  (to  show  the  local  slope  of  the  control  space)  with  the 
number  of  iterations;  (ii)  a  direct  comparison  of  the  control  functions  (the  initial  4>i,  the  target 
4>0  and  the  optimal  4>p).  The  cost  function  value  is  reduced  by  an  order  in  only  2  main  iterations 
(with  several  line-minimization  steps  for  each  main  iteration)  and  by  an  order  of  103  in  12  main 
iterations;  the  local  gradient  norm  is  also  reduced  by  an  order  of  102  and  reaches  a  ‘flat  spot’  in 
12  main  iterations.  At  the  same  time,  the  control  is  changed  from  the  initial  4>i  to  the  optimal  4>p- 
As  we  can  expect  from  the  reduction  of  the  cost  function,  the  optimal  control  matches  very  well 
with  the  target  control  until  approximately  t=8.  The  control  cannot  be  improved  much  after  t=8 
because  of  the  information  delay  due  to  the  distance  between  the  observation  region  fl0,  where 
the  cost  function  is  defined,  and  the  plate,  where  the  control  is  applied.  The  travelling  time 
between  these  two  points  would  require  the  simulation  and  the  cost  function  to  include  future 
events  at  t  >  10  in  ft0  for  optimization  of  the  plate  velocity  in  the  time  period  8<t  <10. 

To  further  confirm  the  accuracy  of  the  gradient  (and  the  associated  efficiency  for  the 
optimization),  we  validate  the  gradient  information  itself  by  a  comparison  between  the  gradient 
obtained  by  the  adjoint  approach  and  the  gradient  computed  directly  by  a  finite  difference 
approach.  The  finite  difference  approach  perturbs  the  control  with  a  small  value  at  each  time 
moment  individually  and  calculates  the  gradient  directly.  Although  the  direct  approach  is 
straightforward  and  considered  accurate,  its  computational  cost  is  proportional  to  the  number  of 
control  parameters.  Figure  36  then  compares  the  gradients  computed  by  the  adjoint  method  with 
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the  ones  computed  by  the  finite  difference  method,  which  shows  good  accuracy  with 
extraordinary  time  saving. 

To  test  the  sensitivity  of  the  algorithm  to  initial  values  and  verify  its  independence  of  a 
particular  function  format,  we  run  the  same  test  with  the  same  target  control  but  with  a  very 
different  initial  control.  As  it  is  shown  in  figure  37,  the  same  convergence  rate  is  observed,  and 
the  same  optimal  solution  is  reached  (for  the  controllable  time  region). 


Figure  35.  Adjoint-based  optimization  for  the  plunging  velocity:  (a)  change  of  the  cost 
function  J  and  the  gradient  norm  Igl  with  the  number  of  iterations;  (b)  comparison  of  the 
initial  control  ( - ),  the  target  control  (j)0  ( —  • — )  and  the  optimal  control  4>p  ( - ). 


Figure  36.  Comparison  of  the  gradient  g  between  the  values  computed  by  the  adjoint  method 
(line)  and  the  finite  difference  approach  through  direct  perturbation  (square). 
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Figure  37.  Adjoint-based  optimization  for  the  plunging  velocity  with  a  different  initial  control. 

3.4  Applications  to  optimize  the  kinematics  and  deformation  of  flapping  wings 

In  this  section,  the  adjoint-based  algorithm  is  applied  to  three  cases:  (i)  a  rigid  plate  in 
combined  plunging  and  pitching  motion  with  an  incoming  flow;  (ii)  a  flexible  plate  in  plunging, 
pitching  and  prescribed  deformation  with  an  incoming  flow;  (iii)  a  three  dimensional  (rigid  or 
flexible)  plate  in  hovering  motion.  For  the  first  two  cases,  the  control  goals  are  drag  reduction 
(i.e.  thrust  enhancement)  and  propulsive  efficiency  improvement;  the  control  parameters  include 
the  phase  delay  between  the  plunging  and  pitching  (for  the  rigid  plate)  and  the  phase  delay  and 
amplitude  of  the  pitching  motion  and  the  first  two  eigenmodes  (for  the  flexible  plate).  For  the 
last  case,  the  control  goal  is  lift  increase;  the  control  parameters  include  the  translational 
amplitude,  the  amplitude  and  phase  delay  of  the  rotational  motion  and  (for  the  flexible  plate)  the 
amplitude  and  phase  delay  of  the  first  eigenmode  in  chordwise  bending. 

1 )  Optimization  of_a  rigid  plunging  and  pitching  plate 

Anderson  et  al.  [61]  suggested  the  critical  role  of  the  phase  delay  between  the  pitching  and 
the  plunging  motion  in  thrust  production  and  propulsive  efficiency,  therefore,  it  is  chosen  to  be 
the  control  parameter  for  drag  reduction.  The  phase  delay  is  optimized  in  two  steps.  First,  we 
consider  the  phase  delay  to  be  a  constant  in  time  and  search  for  its  optimal  value.  Second,  after 
the  optimal  constant  is  found,  we  use  this  value  as  an  initial  condition  to  optimize  the  phase 
delay  further  by  allowing  its  variation  in  time.  The  two-step  strategy  brings  the  control  to  a  good 
neighbourhood  in  a  simple  control  space  (with  a  constant  phase  delay)  before  it  becomes  a  more 
complex  control  space  (with  a  time-varying  phase  delay). 

For  the  constant  phase  delay,  we  start  with  an  arbitrary  initial  value  (p1  =  30°.  It  only  takes  2 
iterations  to  reach  the  optimal  value  cpp  =  —77.3°.  The  corresponding  drag  coefficients,  used  as 
the  cost  function,  is  reduced  from  0.120  to  -0.198.  For  the  time-varying  phase  delay,  the  initial 
value  is  chosen  to  be  the  optimal  constant  phase  delay.  The  optimization  algorithm  reduces  the 
drag  further  by  53%  from  -0.198  to  -0.303  (figure  38). 
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Figure  38.  Optimization  for  the  drag  reduction  of  a  rigid  flapping  plate:  (a)  the  cost  function  for 

drag  reduction  versus  the  number  of  iterations;  (b)  the  optimal  constant  control  ( - )  and  the 

optimal  time-varying  control  ( - ),  with  dash-dot  lines  representing  the  constraints. 


Figure  39  shows  leading-edge  vortex  is  moved  from  the  backside  of  the  plate  to  the  front  side. 
The  new  location  of  leading-edge  vortex  produces  a  low-pressure  region  in  front  of  the  plate 
which  generates  a  large  thrust  force  in  both  the  upstroke  and  the  downstroke. 


Figure  39.  The  vortex  structures  of  a  flapping  plate  with  (a)  the  initial  time  constant,  (b)  optimal 
the  time-constant  and  (c)  the  optimal  time-varying  phase  delay  during  the  upstroke. 

Further  analysis  of  the  optimal  control  indicates  that  the  strategy  of  adjusting  the  phase  delay 
for  thrust  performance  works  to  increase  the  angle  of  attack  magnitude  when  the  plunging  speed 
is  large  to  focus  on  the  enhancement  of  lift-induced  thrust,  and  works  to  decrease  the  angle  of 
attack  magnitude  when  the  plunging  speed  is  small  to  focus  on  the  reduction  of  viscous  drag. 
Meanwhile  the  signs  for  the  angle  of  attack  and  the  plunging  velocity  are  kept  opposite  to  always 
generate  positive  lift-based  thrust. 

When  propulsive  efficiency  is  considered,  the  optimization  result  shows  that  the  reduction  of 
viscous  drag  takes  a  more  dominant  role  because  of  its  impact  on  overall  power  consumption. 
The  mechanism  to  increase  power  output  generally  works  against  the  mechanism  to  reduce  the 
total  power  consumption,  keeping  the  viscous  drag  in  check  becomes  the  key  to  finding  a 
balance  between  the  two  mechanisms  for  optimization. 

2)  Optimization  of  a  flagging  flexible  plate 
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In  this  section,  flexibility  is  added  to  the  flapping  plate,  and  the  analysis  of  its  aerodynamic 
effect  is  undertaken  via  a  comparison  of  the  optimal  and  non-optimal  solutions. 

The  flexible  flapping  plate  has  a  combined  motion  with  plunging,  pitching  and  deformation, 
which  uses  eigenmodes  of  a  cantilevered  Euler-Bernoulli  beam  as  basis  functions.  The  control 
for  a  flexible  flapping  plate  is  then  defined  by 

4*  —  {^0>  <Po.  ^Pn} 

where  a0,  cp0are  the  amplitude  and  phase  delay  of  the  rigid  body  pitching  motion,  and  an,  cpn  (n  > 
0)  are  the  amplitude  and  phase  delay  of  each  eigenmodes  for  deformation  with  n  representing  the 
levels  of  flexibility. 

It  is  shown  in  table  4  and  5  that  the  adjoint-based  optimization  provides  solutions  for  all 
cases  to  change  from  drag  to  thrust  (and  to  larger  thrust).  The  comparison  between  the  rigid  plate 
(case  1)  and  the  flexible  plates  (case  2  and  3)  shows  significant  contribution  from  the  flexibility 
to  the  thrust  improvement.  At  Reynolds  number  300  (Group  I),  the  optimal  flexible  plate  (case  3) 
provides  111%  more  thrust  than  the  rigid  plate  (case  1);  at  lower  Reynolds  number  100  (Group 
II),  the  impact  is  even  more  dramatic  with  1069%  more  thrust  from  the  rigid  plate  (case  1)  to  the 
flexible  plate  (case  3).  For  both  Reynolds  numbers,  a  higher  level  of  flexibility  (case  3,  with  2 
eigenmodes)  provides  more  thrust  than  a  lower  level  of  flexibility  (case  2,  with  1  eigenmode).  It 
is  worth  noting  that  the  thrust  brought  in  by  more  flexibility  comes  with  some  small  cost  of 
efficiency,  since  the  propulsive  efficiency  is  not  part  of  the  cost  function  at  this  moment. 

The  current  study  also  shows  that  flexibility  helps  to  reduce  the  sensitivity  of  the  propulsion 
performance  to  the  Reynolds  number.  For  a  rigid  plate  (n=0),  the  optimal  thrust  is  reduced  by 
86.5%  from  0.26  to  0.035  when  the  Reynolds  number  changes  from  300  to  100;  for  a  flexible 
plate  (n=2),  the  optimal  thrust  is  only  reduced  by  25.5%  from  0.549  to  0.409.  It  suggests  that 
flexibility  helps  to  maintain  the  aerodynamic  performance  in  complex  environment  with  variable 
Reynolds  number  (e.g.  flow  with  gust),  and  flexibility  also  plays  a  more  significant  role  in  lower 
Reynolds  number  region  such  as  flapping- wing  flight  of  insects. 

For  detailed  comparison,  figure  40  shows  the  history  of  drag  forces  for  the  optimal  rigid  plate 
(case  1),  the  optimal  plate  with  small  deformation  (case  2),  the  optimal  plate  with  large 
deformation  (case  3)  and  the  reference  rigid  plate  (case  4,  by  removing  the  flexibility  directly 
from  case  3)  for  Re=300  (Group  I).  Among  the  half-cycle,  the  one  with  small  deformation  (case 
2)  reduces  the  drag  by  a  small  amount  near  t=0  and  t=1.5;  the  one  with  large  deformation  (case  3) 
reduces  the  drag  by  a  large  amount  near  t=0  with  a  slight  increase  near  t=l;  the  reference  (case  4) 
keeps  the  same  or  even  better  drag  reduction  at  t=0  but  cannot  maintain  the  same  saving  later, 
instead  it  adds  large  drag  for  let  <2. 

Figure  41  shows  visually  the  flapping  kinematics,  deformation  and  vortex  structures  for 
better  understanding  of  the  control  mechanism.  Flexibility,  especially  large  deformation,  allows 
the  plate  to  have  a  large  angle  of  attack  and  hold  the  leading-edge  vortex  to  the  front  for  the 
benefit  of  lift-induced  thrust,  while  the  overall  profile  stays  small  (by  deformation)  to  avoid  an 
increase  in  viscous  drag.  It  is  obvious  that  the  reference  case  4,  which  keeps  the  same  pitching 
(i.e.  angle  of  attack)  at  the  leading  edge  but  removes  all  the  deformations,  has  the  same  benefit 
from  the  lift-induced  thrust  (near  t=0)  but  adds  a  huge  amount  of  viscous  drag  later. 

Table  4.  Optimization  results  for  drag  reduction  with  different  flexibility  at  Re=300:  cases  0  has 
the  initial  control;  cases  1,  2  and  3  have  the  optimal  controls  with  different  levels  of  flexibility; 
case  4  is  a  reference  case,  which  removes  all  flexibility  terms  from  case  3  and  keeps  the  rigid 
plunging  and  pitching  motion  the  same.  r|  is  the  propulsion  efficiency. 
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Group 

Cases 

n 

cd 

i] 

0 

0 

0.138 

-0.0373 

1 

0 

-0.260 

0.184 

I 

2 

1 

-0.376 

0.296 

3 

2 

-0.549 

0.273 

4 

0 

0.696 

-0.141 

Table  5.  Optimization  results  for  drag  reduction  with  different  flexibility  at  Re=100. 


Group 

Cases 

n 

Cd 

l] 

0 

0 

0.165 

-0.473 

1 

0 

-0.035 

0.0194 

II 

2 

1 

-0.177 

0.0892 

3 

2 

-0.409 

0.116 

Figure  40.  The  history  of  drag  forces  in  one  period  for  case  1  ( - ),  case  2  ( - ),  case  3  ( —  • 

— )  and  case  4  ( )  of  Group  I. 
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Figure  41.  Comparison  of  kinematics,  deformation  and  vortex  structures  of  (a)  the  rigid  plate 
(case  1),  (b)  the  plate  with  small  deformation  (case  2),  (c)  the  plate  with  large  deformation  (case 
3)  and  (d)  the  reference  (case  4),  with  Re=300  (Group  I),  at  t=0,  1  and  1.5. 


Adjoint-based  optimization  also  improves  the  propulsive  efficiency  substantially  (table  6). 
The  optimal  solution  shows  that  plate  tries  to  blend  into  a  streamline  shape  to  reduce  the  viscous 
drag.  Meanwhile  holding  a  large  leading-edge  vortex  is  no  longer  an  option  with  its  cost  on 
viscous  drag.  Same  as  that  in  drag  reduction,  flexibility  also  largely  reduces  the  sensitivity  of 
propulsive  efficiency  to  Reynolds  number  and  allows  a  flapping  wing  at  low  Reynolds  number  to 
enjoy  the  high  aerodynamic  performance  existing  in  the  higher  Reynolds  number  regime. 


Table  6.  Optimization  results  for  propulsive  efficiency  with  different  flexibility  at  Re=300. 


Group 

Cases 

n 

Ci 

V 

0 

0 

0.138 

-0.0373 

1 

0 

-0.223 

0.247 

III 

2 

1 

-0.307 

0.373 

3 

2 

-0.347 

0.377 

3}  Optimization  of_a  three-dimensional  plate  in  hovering  motion 
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In  this  section,  we  extend  the  study  from  a  two-dimensional  plate  to  a  three-dimensional 
plate.  We  choose  the  same  rigid  plate  in  hovering  motion,  which  has  been  studied  extensively  by 
Trizila  et  al.  [44],  as  our  base  case.  With  initial  numerical  simulations  combined  with  surrogate 
modelling,  Trizila  et  al.  [44]  achieved  the  complete  map  of  a  three-parameter  design  space 
including:  the  amplitude  of  translational  motion  x0,  the  amplitude  0a  and  the  phase  delay  of 
rotational  motion  (p.  Such  a  space  exploration  is  important  for  some  cases.  However,  an  adjoint- 
based  optimization  could  be  much  more  efficient  if  optimization  is  the  major  concern  and  the 
result  is  more  accurate  without  involving  a  surrogate  model. 

Table  7  shows  respectively  the  initial  control  and  its  optimal  value  after  4  main  iterations. 
The  lift  coefficient  is  improved  from  0.338  to  0.620  by  83.4  %.  A  closer  comparison  shows  that 
our  optimal  values  match  well  with  the  high-lift  region  shown  in  the  contour  map  provided  by 
the  design  space  exploration  through  surrogate  modelling  in  Trizila  et  al.  [44]. 

Table  7.  Comparison  of  the  initial  and  the  optimal  controls  and  the  resulting  lift  coefficients  for  a 
rigid  hovering  plate. _ 


Control 

x0 

0« 

c, 

Initial 

1.5 

O 

O 

v£> 

90° 

0.338 

Optimal 

1.0 

45° 

120° 

0.620 

Figures  42  compares  the  time  history  of  lift  coefficients  with  the  initial  and  the  optimal 
controls.  Compared  to  the  initial  control,  the  plate  with  the  optimal  control  has  a  larger  angle  of 
attack  at  t=0.  This  generates  a  stronger  leading-edge  vortex  (figure  43),  which  causes  large 
pressure  difference  in  upper  and  lower  surface  of  the  plate  (figure  44)  and  produces  larger  lift.  At 
t=0.2,  in  the  case  with  the  initial  control,  the  vortex  pair  generates  a  momentum  on  the  plate’s 
upper  surface  (figure  43),  therefore  it  produces  downward  force  (figure  44).  On  the  other  hand, 
the  plate  with  the  optimal  control  has  an  advanced  rotation.  The  momentum  induced  by  the 
vortex  pair  acts  on  the  plate’s  lower  surface,  thus  generating  lift  instead.  These  results  suggest 
that,  in  hovering  flight,  both  the  plunging/pitching  and  the  wake  capture  are  important  for  lift 
generation.  The  timing  of  the  plate  interacting  with  the  wake  is  closely  related  to  the  phase  delay. 
The  vortex  structures  depicted  in  figure  45  by  an  isosurface  of  Q  criterion  clearly  show  much 
stronger  downwash  from  the  optimal  motion  for  the  increase  of  lift  performance. 

When  chord-wise  deformation  is  considered  during  the  optimization,  the  lift  coefficient  is 
pushed  from  0.62  to  0.88  for  another  42%  increase.  The  flexibility  is  able  to  push  the  lift 
coefficient  higher  by  ‘cupping’  more  fluid  on  the  lower  surface  and  deforming  to  larger  vertical 
projection  profile  which  lead  to  an  even  stronger  downwash  as  well  as  larger  pressure  difference 
between  the  upper  and  lower  surfaces  of  the  plate,  both  for  lift  benefit. 
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Figure  42.  The  time  history  of  the  lift  coefficients  for  a  rigid  hovering  plate  with  the  initial 
control  ( - )  and  the  optimal  control  ( - ). 


Figure  43.  The  z-component  of  vorticity  at  the  middle  plane  for  a  rigid  hovering  plate  with  (a) 
the  initial  control  and  (b)  the  optimal  control  at  different  time  moments. 
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t  =  —0.2  ?  =  0  f  =  0.2 


Figure  44.  The  pressure  distribution  of  a  rigid  hovering  plate  with  (a)  the  initial  control  and  (b) 
the  optimized  control  on  the  lower  (left)  and  upper  (right)  surfaces  at  different  time  moments. 


/  =  —0.2  ?  =  0  f  =  0.2 


Figure  45.  The  vortex  structures  of  a  rigid  hovering  plate  with  (a)  the  initial  control  and  (b)  the 
optimal  control  at  different  time  moments. 
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